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Abstract. We review recent progress on constructing non-equilibrium steady state 
density operators of boundary driven locally interacting quantum chains, where driving 
is implemented via Markovian dissipation channels attached to the chain’s ends. We 
discuss explicit solutions in three different classes of quantum chains, specifically, the 
paradigmatic (anisotropic) Heisenberg spin-1/2 chain, the Fermi-Hubbard chain, and 
the Lai-Sutherland spin-1 chain, and discuss universal concepts which characterize 
these solutions, such as matrix product ansatz and a more structured walking graph 
state ansatz. The central theme is the connection between the matrix product form of 
nonequilibrium states and the integrability structures of the bulk Hamiltonian, such 
as the Lax operators and the Yang-Baxter equation. 

However, there is a remarkable distinction with respect to the conventional quantum 
inverse scattering method, namely addressing nonequilibrium steady state density 
operators requires non-unitary irreducible representations of Yang-Baxter algebra 
which are typically of infinite dimensionality. Such constructions result in non- 
Hermitian, and often also non-diagonalisable families of commuting transfer operators 
which in turn result in novel conservation laws of the integrable bulk Hamiltonians. 
For example, in the case of anisotropic Heisenberg model, quasi-local conserved 
operators which are odd under spin reversal (or spin flip) can be constructed, whereas 
the conserved operators stemming from orthodox Hermitian transfer operators (via 
logarithmic differentiation) are all even under spin reversal. 
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1. Introduction 

1.1. The key concepts of quantum integrability 

Exact solutions of nontrivial yet simple physical models are of paramount importance 
in statistical and quantum physics. On one hand, an exactly solved model often 
characterises universal behaviour of a more general class of possibly unsolvable models 
and thus represents the best possible exact understanding of physical reality. On the 
other hand, the tricks developed in the course of deriving such solutions have often lead 
to general development in mathematical methodology or even to novel mathematical 
concepts. Famous examples being for instance, Hans Bethe’s solution of the Heisenberg 
model of magnetism in 1931 and Lars Onsager’s solution of planar (2D) Ising model 
in 1944. Remarkably, these two threads have merged in the works of C. N. Yang [95] 
and Rodney Baxter [5], giving birth to the celebrated Yang-Baxter equation (YBE), 
also known as the star-triangle equation, representing the most general characterisation 
of integrability known to date. Moreover, abstract algebraic characterisation of YBE 
lead Vladimir Drinfeld in 1980’s to introduce the concept of quantum groups and 
develop their representation theory, together with Jirnbo, Reshetikhin, Sklyanin and 
many others. 

An important ingredient of the theory of quantum integrability is the concept of 
auxiliary Hilbert space which interacts with the physical quantum Hilbert space via 
the so-called Lax (or scattering) operator, solving the YBE together with the so-called 
R-matrix which represents the ‘internal’ scattering between a pair of auxiliary spaces. 
Considering explicit matrix representations pertaining to finite dimensional auxiliary 
space — typically being a fundamental representation of the corresponding quantum- 
or Lie-group symmetry — resulted in the very successful quantum inverse scattering 
method [S7J:2S] for diagonalising integrablc many-body Hamiltonians and computing 
their equilibrium correlation functions, also known as the algebraic Bethe ansatz 
(ABA). Remarkably, literally the same technique has beed adapted to solve certain 
nonequilibrium classical driven diffusive many-body systems in one-dimension, namely 
the so-called simple exclusion processes (SEP), for a simple reason that their Markov 
chain generator can be identified with the Heisenberg-like Hamiltonian. However, 
ABA calculations often result in an implicit solution written in terms of a coupled 
set of algebraic (or transcendental) equations [the Bethe equations (BE)], or in the 
thermodynamic limit (TL), in terms of coupled integral or functional equations (so- 
called Bethe-Takahashi equations). 

In a seminal paper [19], Derrida and coworkers have been able to circumvent this 
problem by writing an explicit solution for the steady state of symmetric and asymmetric 
simple exclusion processes (ASEP/SSEP) in terms of a matrix product ansatz (MPA). 
MPA in turn also allowed for explicit calculation of all physical observables and 
correlation functions in the nonequilibrium steady state (NESS), see e.g. [T8l [821 18] for 
review. It is notable however, that a particular MPA appeared earlier as a ground state of 
a valence bond solid in one-dimension (AKLT model p]), or as a convenient family of the 
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so-called finitely correlated states of quantum spin chains ra. and nowadays represents 
a cornerstone of the density-matrix-renormalization-group (DMRG) method, a state-of- 
the art technique for simulation of strongly correlated systems in one-dimension |81j . It 
has later also been recognised that ABA method for equilibrium integrable Heisenberg 
spin chains can be re-phrased in terms of MPA for all eigenstates, from which BE can 
be equivalently derived mmm- 

However, the problem of integrability in the combined paradigm of driven diffusive 
systems and quantum many-body interactions resisted until 2011 when NESS density 
operator of the boundary driven Lindblad master equation for the anisotropic Heisenberg 
spin 1/2 chain (XXZ model) has been solved in terms of MPA [69) [TO] by the author 
of the present topical review]]]] The solution appeared disconnected from conventional 
theory of quantum integrability at the first sight, but has later [771EZIE1] been related 
to infinite-dimensional solutions of YBE pertaining to non-unitary (or in mathematical 
terminology, non-integrable ) irreducible representations of quantum group symmetry. 


1.2. The purpose and summary of the review 


We have by now managed to derive a number of exact MPA solutions of NESS of 
boundary driven quantum chains for different types of integrable locally interacting 
bulk Hamiltonians and different boundary dissipators (diffusive driving). The unifying 
picture of such non-equilbrium quantum integrability has been emerging slowly from 
studying various quite specihc situations scattered over several rather technical papers, 
so it is perhaps now a good moment to wrap these results within a single topical review 
article. The purpose of the present paper is thus to present a coherent and up-to-date 
review of the progress on solving the challenging problem of integrable boundary driven 
diffusive (or better to say, disipatively driven) quantum systems. The problem may be 
of quite general interest for mathematical physics and statistical physics community as 
it represents in some sense a minimal description of dissipative or incoherent driving of 
a quantum many-body system without affecting coherent macroscopic character of its 
bulk. One can view this approach as an incoherent or dissipative analogue to a standard 
trick in nonequilibrium transport of coherently driven systems, where one introduces a 
bias in electro-chemical potential replacing a real electric held in the bulk. In contrast 
to most of literature on integrable systems, this review takes what might be considered 
a bottom-up approach. Namely, we first demonstrate how to explicitly construct various 
exact physical solutions and only later think or elaborate on their abstract mathematical 
properties and meaning. 


In the following subsection |1.3| we physically motivate the paradigm of boundary 
driven quantum master equation of the Lindblad form. In section [2] we then describe 
the MPA solution of NESS in case of the simplest generic integrable model, specifically 


f The term ‘driven diffusive systems’ can in fact be misleading in our context since, as we shall see 
later, the competition between dissipation and coherent quantum many-body interactions can lead to 
a plethora of transport behaviours, ranging from ballistic, anomalous, diffusive, to insulating. 
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the XXX or more generally, the XXZ spin 1/2 chain. Having the explicit MPA form, 
we also discuss in detail explicit analytic computation of observables and correlation 
functions in NESS, and its relation to the quantum group U q (sl 2 ) and the corresponding 
solutions of YBE. In sections 0]and[5]we then discuss MPA solutions in two other families 
of boundary driven systems, specifically in the one-dimensional fermionic Hubbard 
model and in Lai-Sutherland spin-1 chain, respectively. The former is intriguing as 
the corresponding Lax operator generating the MPA seems not to be connected to 
obvious symmetries of the model neither to the celebrated Shastry’s i?-matrix, while 
the latter is fundamentally interesting since the NESS is macroscopically degenerate and 
can be parametrised by a thermodynamic variable yielding a nonequilibrium analog of 
the grand-canonical ensemble. In section [6] we discuss the relevance of our exact far- 
from-equilibrium results for linear response physics of quantum transport, such as the 
existence of novel quasilocal conserved operators and rigorous lower bounds on zero 
frequency Green-Kubo transport coefficients. In section [7] we outline a subset of most 
exciting and urgent open problems and conclude. 

We stress that the present review contains also several original results w.r.t. the 
existing literature. In particular, a mixture of asymmetric coherent (boundary fields) 
and incoherent (boundary dissipation) has been worked out completely for XXZ chain 
in subsection 2.6.1, as well as a fully analytic calculation of the nonequilibrium partition 
functions in the isotropic XXX chain for a variety of asymmetrized boundary drivings 
in subsection 13.21 


1.3. The paradigm of boundary driven quantum master equation 


Let us consider a finite d— dimensional local physical Hilbert space TL P — C d , e.g. of a 
quantum spin s = (d — l)/2, and consider a quantum chain of length n defined on a 
tensor product space TL® n = (££)™ =1 with a Hamiltonian which can be written in terms 
of a sum of local interactions h XtX+ i acting nontrivially only over a pair of neighbouring 
sites (x, x + 1) 


72—1 

H ^ ^ hxjX+l* 

X=1 


( 1 . 1 ) 


Within the theory of open quantum systems mi (see Refs. [3, [16] for more 
mathematical accounts) incoherent markovian quantum dissipation can be completely 
described by a set of quantum jump operators { L M E End("Hp n ); p = 1, 2 ...}, also called 
Lindblad operators, so that the evolution of the system’s many-body density operator 
pit) satisfies the Lindblad-Gorini-Kossakowski-Sudarshan master equation [51, 20] 

^p(t) = Cp(t) := -i [H, p(t )] + (2 L nP(t)Ll - {L\L^ p(t)}) . (1.2) 


Eq. (1.2) defines a general family of markovian dynamical semigroups lA{t) = exp(LC), 
U(t)U(t') = lift + £'), t,t' E M + , which preserve hermiticity, positivity and trace of the 
density matrix p{t) at all times. In fact, the semigroup (1.2) can be derived from the 
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Figure 1. Illustration of a dissipatively boundary driven quantum spin chain where 
the first and the last particle of the chain are immersed into a pair of distinct markovian 
quantum baths corresponding to different values of thermodynamic potentials. 


microscopic unitary evolution of the universe = system <g) environment, provided the 
following conditions are met: (i) the initial state of the universe is factorized p(0) <g)p e nv, 
(ii) the coupling between the system and the environment is weak so the second- 
order Born-Dyson series can be used, and (iii) the dynamical correlation functions of 
environment observables that are coupled to the system decay on time scales for which 
the system’s dynamics can be considered as frozen ( secular approximation). 

We shall furthermore assume that the incoherent quantum processes are ultra-local, 
meaning that all jump operators are of the form 


L,= 


1 (Jri — 1 


or 


L/J, — 1 d™- 1 




4 G End('H p ), (1.3) 

i.e., they act nontrivially either on the left or the right boundary of the chain (Fig |T|). 
Within the microscopic derivation mi , this additional assumption is generically justified 
only if the on-site part of h XjX+ \ is much larger (in operator norm) than the genuine 
interaction part [H], meaning that a local disturbance on the boundary site does not 
spread to the interior before it gets dissipated so the dissipator can be assumed to act 
locally. However, there is an alternative phenomenological description of the boundary 
driven quantum master equation ( 1.2|1.3 ) in terms of the so-called repeated interactions 
protocol 3j). CE] which is free from any assumptions about the microscopic dynamics. 
In this protocol it is assumed that end spins/particles of the chain are put in contact 
(interaction) for a short amount of time 5t with a pair of auxiliary spins/particles which 
are assumed to be prepared in two different thermal, canonical, grandcanonical, or any 
other equilibrium states. In the next time interval 5t, the auxiliary spins, the states 
of which may have already slightly changed, are replaced by a fresh, independent pair 
of thermal auxiliary spins, and so on. One may imagine two running belts carrying 
thermally prepared auxiliary spins in separable states and moving along each side of 
the chain at some speed while interacting with the ends of the central chain. Writing 
dynamical evolution of the central chain density matrix in the high-speed limit 5t —>• 0 
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one precisely recovers (1.2). It should be noted that closely related protocols could 
nowadays be implemented using contemporary cold atom experiments. 

We stress, however, that our main motivation for studying the problem ( 1.2|1.3 ) is 
in its conceptual simplicity and mathematical elegance of formulation as the minimal 
model which can describe nonequilibrium driving of a coherent many-body system in 
one dimension. Such a setup allows for an efficient numerical (DMRG) simulation of 
quantum transport for generic non-integrable bulk interactions [78] [SB]. Even though 
this review focuses on the case where the bulk Hamiltonian is strongly interacting, 
we start by mentioning that in the non-interacting case, where the Hamiltonian is 
usually mapped to XY spin 1/2 chain and the Lindblad jump operators are linear in the 
corresponding Wigner-Jordan fcrmi operators, one is able to calculate analytically the 
full NESS density operator, all its observables, as well as the relaxation dynamics of the 
master equation (1.2), in terms of nonequilibrium dissipative quasi-particle excitations 


Furthermore, one is able to write exact solutions in the non-interacting case 
even in some situations where the jump operators are quadratic but Hermitian, such as 
for the so-called dephasing noise, which for exact solvability has to be homogeneously 
distributed in the bulk of the chain |J5j. This solution can in fact be written in terms 
of a simple MPA El with 4x4 auxiliary matrices and allows for some further solvable 
generalisations [23]- Such models can be further generalized as hybrid quantum-classical 
rnarkov chains El, where the incoherent part of dynamics exactly coincides with the 
classical SEP. 

Focusing on the concept boundary of driven many-body systems in this review, 
where the baths are described only effectively, we have to refrain from discussing a bulk 
of related and also highly topical literature on nonequilibrium steady states with infinite 
(microscopically formulated) baths. 


2. Anisotropic Heisenberg spin-1/2 chain 

We shall start by considering arguably the simplest integrable model with strong 
interactions — the XXZ model, namely a homogeneous spin 1/2 chain (d = 2) with 
nearest neighbour anisotropic Heisenberg interaction with anisotropy parameter A, 

h = 2cr + 8) a~ + 2cr _ <8 a + + Acr z 8 cr z , (2.1) 

where a ± = |(<r x ± icr y ) and <7 x ,cr y ,cx z are the usual 2x2 Pauli matrices, so that the 
Hamiltonian ( |1.1[ ) can be written as 

n—1 

H = ^ l 2 x-i 8 h 8 ton-x-i. (2.2) 

X=1 

Embedding the Pauli operators into End("Hp n ), cr" := l 2 *-i 8cr“8 translationally 

shifted interactions read h XjX+1 = 2(J Z°x+l + 2a x °t+i + A(T X+i- 

We start by considering the simplest nontrivial dissipative driving with only two 
jump processes coupled to bulk unitary dynamics at a single dissipation rate £ E M + , 

L\ = , Lo = 'f£cr n . (2.3) 
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As the total z-component of magnetization M = V^"_ 1 o x and hence the number of 
up-spins are conserved in the bulk 

[H, M] — 0, (2.4) 


the incoherent processes (2.3) can be interpreted as a pure source of up-spins on 


the left end and a pure sink of up-spins on the right end. NESS density operator 
Poo = lirn^oo exp(t£)p(0) can be considered as a fixed point of the propagator, or null 
eigenvector of the Liouvillian 

E Poo = 0, where C = — i&dH + £'D (T + + e'D^-, (2.5) 

where the Lie derivative map (ad H)p := [H, p] and an elementary dissipator map 

V L {p) = 2LpL t -{L%p} (2.6) 

have been introduced. 


2.1. Uniqueness of NESS 


Let us first show EB that under quite general conditions, the Liouvillian ( |2.5[ ) possesses 
a unique NESS, i.e. the fixed point is independent of the initial state p(0). 

We start by noting a theorem due Evans 129 and Frigeiro [28] (with the subject 
nicely reviewed by Spohn [HS] ) which essentially states that NESS is unique iff the set of 
operators {H, Li, L\,L 2 , L\ ■..} generates, under multiplication and addition, the entire 
algebra End("Hp n ) of operators over a quantum chain on n sites. Indeed, this is easy 
to demonstrate explicitely even considering only a triple of operators {H, erf , erf} while 


uniqueness then trivially extends to all cases including, and generalizing, (2.3) where 
the set {L^, Ljp p = 1, 2 ...} contains, up to scalar prefactors, either a pair of, of, or a 
pair of, of due left-right inversion symmetry. 

Namely, one observes the following recursive operator identities: 


<*t = i H ^ i]]> 


o' = 


a 


x—2 




x — 3,4 ..., n, 


(2.7) 

( 2 . 8 ) 


which generate the entire set {erf ; x — 1,..., n} starting from just H and of. Similarly, 
{er~;x = 1,n) are generated by Hermitian adjoints of Eqs. ( 2.7[2.8 ) starting from 
H and of. The set {of, of]x = 1,..., n} then generates all elements of End('H® n ) by 
multiplication and addition. 

For related recent general results on (non)uniqueness of fixed points and 
characterization of the space of steady states of Lindbladian dynamics, see, e.g., 

Refs. [41 SS]. 


2.2. Matrix product solution - isolating defect operator method 

We shall now present an ad hoc method which generates the MPA of NESS fixed point 
Poo for the XXZ model following Ref. [70] (also [69]), which we term an isolated defect 
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operator (IDO) approach. Later we shall in subsect. 2.4 rederive this solution in a more 


elegant way using a non-Hermitian infinite-dimensional Lax matrix of the XXZ model, 
but we believe it may be of interest to investigate both approaches. For instance, if for 
some model the appropriate Lax matrix for the problem is not known, IDO strategy 
may sometimes be used to devise a systematic method for determining MPA by an 
automated procedure (see e.g. Ref. can- 

Let us first show that NESS density operator allows a particular factorization in 
terms of a non-Hermitian amplitude operator D n (e). 

Lemma 1. Let G End('Hp"') satisfy the following conditions (defining relations): 

(i) a recursion identity for the bulk, setting hR := a 0 := t 2 , 

[. H , Q n ] = -E(a z <g O n _i - O n _i <g a z ) 

and (ii) the recursion identity for the boundaries 

fl n = a 0 g n„._i + a + g 

for some unspecified operators e End(%p ^ n_1 ' ) ) 

Roc 


(2,9) 


Poo 


tr R r 


Roc 


a u + D n _ 1 0 a (2-10) 

Then, the density operator 

( 2 . 11 ) 


satisfies the fixed point condition (2.5) 


Proof. We need to prove that £(D n Di) = 0, i.e. 


1£ 


-1 


[H, n»nt] = v „+ (nX) + A- (nX)- 


( 2 , 12 ) 


Using the Leibniz rule for Lie derivative, and (i), LHS of ( |2.12 ) can be transformed to 

{i£~ l [H, n n ]) f + n n ]nt = 

D n (cr z g> ST2^_i) - D n (D|(_ 1 <g) CT Z ) + (cr z <g) D n _i)Dl - (D n _i g> cr z )Dl. (2.13) 

Applying, respectively, the first and the second identity of (ii) on the first two and the 


last two occurrences of Ll n in (2.13), we obtain 

2a z <g - a + <g 


cr 


O n _ 1 (D+_ 1 ) 1 


20 n _ 1 D^_ 1 g a z - D n _ 1 D),_ 1 (g a - D n _ 1 (D n _ 1 ) t 


a 


(2.14) 


The first term on the RHS of (2.12) can again be transformed by applying the first 
identity of (ii) to 

V a+ (a°) <g SK-itol-i + Va+{°~) ® D n _ 1 (n+_i) t + 

V a+ {a + ) g + A + (cr + a-) g n+_ 1 (n+_ 1 ) t , 

which, by observing the complete local action of the dissipators 

V^(a°)=2a', ■D„,(o ± ) = -o ± , V„ + (o + a~) = 0, 

A-(0°) = -2V, A-(°*) = ~0 ± , A-(°~°*) = 0, 


(2.15) 

(2.16) 
(2.17) 


result in exactly the first three terms of expression (2.14). Analogously, the second term 
on the RHS of (2.12) results in the last three terms of (2.14). □ 
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Note that the condition (ii), Eq. ( 2. 1Q| ) , implies that Q n is unit diagonal and upper 
triangular matrix in the joint eigenbasis {\v \,..., v n ) ; u x G {0,1}} of {cr^}, <J X \p) = 
(1 — 2v x ) |u), where the ordering is determined by the binary ord(V) := Y^=i Pz^ -1 , 

(n| Q n | u) = 1 , (u\ Q n |z/) = 0 if ord(z/) > ord(z/'), (2.18) 


which can be shown by a straightforward induction. The factorisation of NESS (2.11) 
can thus be considered as a (reverse) many-body Cholesky decomposition. In the 
canonical Cholesky decomposition, though, the matrix of f2 n would have to be lower- 
triangular, but this can be trivially mended by considering a spin-reversed problem. 
Namely, noting the spin-reversal symmetry of the Hamiltonian 

0 1 ' 


PHP - 1 = H, P = P 1 = Y[ 


X=1 


1 0 


(2.19) 


the spin-reversed density operator 

R' 0C = PR X P- I = n' n ((l'j, (l' n = Pfl n p-\ ( 2 . 20 ) 

where the matrix of Vt' n is unit diagonal and lower triangular , solves the reverse NESS 
fixed point condition with the source and sink being swapped: 


L' = PLiP” 1 = 


J l ~ 1 ^1 1 ~ V £<J 1 j ±J 2 

We shall now define an abstract auxiliary Hilbert space "H a , with a particular state 
|0) G "Ha, and postulate an MPA for the amplitude operator [§] 

tt n = J2 (0| A ai A a2 • ■ • A an |0) a ai ® <r “ 2 <£>••• a a ", (2.22) 

with a yet-to-be specified tripple of matrices A±, A 0 G End(‘H a ). Throughout this paper 
we write in roman-bold letters operators which act non-trivially, i.e. not as scalars, over 
the auxiliary space P a . We should note that the terms with cr^ have been omitted in 
the ansatz (2.22) which is the key to the solution of the problem. Let us call cr z a defect 
operator here. Looking at the pair of sites without a defect, one finds that commutation 
with Hamliltonian density produces exactly one defect, either on the left or on the right 
tensor factor 

a,a' ft 


L' = PL^P - 1 = 


ea ‘ 


( 2 . 21 ) 


[h, o a ® a a ‘] - ^ 

^s{+,—,0} 


V 


cr 


z , ,Q 7. 

+ 7/3 ° 


o 


(2.23) 


a, of G {+, —, 0}, with the structure constants 7 ^’ a , writing out only the non-vanishing 
elements: 


7 ±’° = ±2A, 


o,± 

7± — "F2, 


7o ±,=F = ±1- 


(2.24) 


The commutator [ H , fl n ] and the entire dehning relation (2.9) should then contain Pauli 
operators with exactly one defect a x . Considering all the terms where the defect operator 


§ An impatient reader should here be directed straight to subsect. 


2.3 
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appears in the bulk 1 < x < n, a sufficient condition for cancelation involves a projection 
on a triple of sites (x — 1 , x, x + 1 ) 

tr ^cr _/3 (8) cr z 8 a~^'[h 8 ) 1 2 + I 2 8 h, cr a <g) cr"' 8 cr""] j A„A a iA a « = 0, (2.25) 

a,a / ,a' , E{+,—,0} 

or, equivalently, in terms of the structure constants 7 ^’“ 

Y {frp’ a 'A a A a/ A p , + ^',’ a AgA a A a ^ = 0 , (2.26) 

a,o'e{+,—,0} 

which represent 9 homogeneous cubic identities for the matrices A±,Aq for fd,f3' G 
{+, —, 0}. Writing these explicitly, we find that only 8 of them are linearly independent: 


[A 0 , A±A t ] = 0, {A 0 , A|} = 2 AA ± A 0 A±, 2 A[A^, A±] = [A T , A|], 

2 A{Ag, A±} - 4A 0 A ± A 0 = {A T , A|} - 2A±A T A±. (2.27) 


Considering the remaining terms of the defining relation (2.9) where the defect operator 
sits at the boundary we obtain two triples of sufficient conditions projecting (in analogy 

ir of sites near 1 

P £ {+) — , 0 } 

(2.28) 


to (2.25)) on cr z 8 cr“ for x — 1, or cr" 8 cr z for x = n, for a pair of sites near the boundary 

Y 7 ^’ a ( 0 | A q A q / = — ier < 0 | Ay, 

a,a'e{+,—, 0 } 

Y r yp’ a A « A «' 1 °) = - ' i£ Ap |o). 


a,a'£{+,—,0} 


In order to fulfil the second defining relation (2.10) with the MPA (2.22), simple 
additional sufficient conditions are 


(0| A_ = 0, A + |0) = 0, (0|A 0 = (0|, A 0 |0) = |0). (2.29) 

We have thus shown that a representation of A±, A 0 fulfilling the cubic bulk relations 
(2.27) together with quadratic-linear boundary conditions ( 2.28|2.29 ) will provide MPA 
for NESS of XXZ chain. Considering |0) as a highest weight state , we define a tower of 
states | k) spanning an infinite-dimensional auxiliary representation space 

H a = lsp{|fc) := (A_) fc |0), k = 0,1, 2 ...}. (2.30) 

Using Dirac’s notation with the dual basis (k |, satisfying (k\k') = 8k,k'i the auxiliary 
operators then naturally take the following tridiagonal form 

OO OO OO 

A 0 = y> fc \k)(k\, A + = Y bk \ k )( k + 1|, A ~ = Y; \ k + ’ ( 2 - 31 ) 


k =0 


k =0 


k =0 


where the 8 relations (2.27) become equivalent to a pair of recurrence relations for 07 , hk 

(2.32) 

(2.33) 


afe+i 2Aafc ~\~ ak —1 0, 

bk+i ~ bk = 2a/c_(_i (Acifc+i — ak) 

with boundary conditions ( 2.28|2.29 ) yielding the initial conditions for the 
a 0 — 1, ai = A + b 0 = ie. 


recurrence 


(2.34) 
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a 

0 


a 

-0 = 

1 


= 8 = 

2 


a 

-O^ 

3 


o 

-0 = 

4 


Figure 2. A semi-infinite graph showing the allowed transitions for building up the 
matrix product operator f2 n (2.371 for the XXZ chain. Brown, red, blue edges e 
denote the index function values w(e) = a°,a + ,a~, respectively. 


The solution can be expressed in terms of Chebyshev polynomials in A, or more 
compactly, re-writing the anisotropy parameter 

cost] := A, (2.35) 


as 


ak = cos(kr]) + 


ic sin(fcr/) 


b k = sm((k + 1)77) 


sin 77 

ic 


Sill 7] 


cos (rjk) 


1 + 


4 sin 2 r) 


sin (rjk) 


(2.36) 


MPA (2.22) can be given another appealing interpretation. Associating the 


auxiliary states | k) with vertices of a directed graph, V = {0,1, 2,...}, and defining a set 
of edges encoding all possible transitions £ = {(k,k), (. k , k+ 1), (k+ 1, k); k — 0,1,...} the 
amplitudes (0| A ai A a2 ■ ■ • A an |0) can be written in terms of a collection of products of 
transition amplitudes along all possible n-step walks , i.e. sequences of connecting edges 
starting at vertex 0 and ending back at 0 in exactly n steps, W n (0, 0) = {ei, e 2 ..., e n G 
£',p{e 1) = 0 ,q{ej) = p(e j+1 ),q(e n ) = 0}, where e = (p(e),g(e)), namely 


f -n 


(ei,...,e„)ew;„(0,0) 


*e n 


cu(ei) <8) ai(e 2 ) <8) • • -u;(e n ). 


(2.37) 


The amplitudes are encoded as a^,k) = a k, a(fc,fc+i) = b k , a^+i,k) = 1, and cu : £ —>• 
End('Hp) is what we shall call an index function which associates a local physical operator 
with each edge of the graph. In our case u(k, k ) = cr°, cj(k, k + 1) = cr + , cu(k + l, k) = cr~, 
where the defect <r z is not in the image of u (see Fig. [2]). Expression (2.37) shall be named 
a walking graph state (WGS) representation, and is a useful concept encapsulating the 
locality of infinite-dimensional MPA. 


2.3. Lax operator for a complex q— deformed spin and some notation 

So far we followed a pedestrian approach and have not used the powerful quantum group 
U q (sl 2) symmetry [42j of the XXZ spin chain [60J which is deeply rooted behind the 
integrability of all equilibrium problems for the model. A fundamental characterisation 
of this symmetry can be written (see e.g. Refs. (25, 145] ) in terms of the so-called RLL 
relation, a version of the YBE over a tensor product triple "H a ® T-L p ® Lip 

Ripipi — <^2)1^1)14^) = Li((^2)L2(<^i)-Ri i2 (<^i — 952)- ( 2 . 38 ) 

Remember that bold-roman letters denote symbols acting (nontrivially) over auxiliary 
space "H a while indices denote the label of the physical space. Here R(<p) is, up to a 
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permutation, the standard trigonometric 4x4 six-vertex R-matrix which we choose to 

write in terms of spin operators over T~L P ® Up 

. sin cp., . l + cosw 1 — cos 03 , . 

R{<p>) — (h+t cos rj) --—— Ism 77 -)-- ——a sm? 7 , (2.39) 


T-Lp) has a universal U q (sL) 


a 


2 " 2 

while the Lax operator (or L-operator) L G End("H 
symmetric form: 

t (, n ( sin (^ + (sinr/)s; \ _ v- a( 

’ V ( sin r l) s f sin ((p-nal)) ££ 

where J = {+, — , 0, z}, and L“ G End('H a ) are its physical space components 
L°(</?, s') = sin ip cos(r]S z s ), L z (<p,s) = cos ip sm(r]S z s ), L ± (<p,s) = (sin77)s;f 
For mathematical applications of quantum group symmetry to abstract construction 
of L-operators and other quantum integrability concepts, see e.g. Refs. PHD]. The 


(2.40) 


(2.41) 


RLL relation (2.38) is in fact equivalent to a complete set of commutation relations 


for the generators s(f, of a q —deformed angular momentum algebra with deformation 
parameter q = e lv 


~ [2 S s ]q — 


sin(277s^ 
sin 77 

,-i\ 


K, sf] 


= ±s^ 


(2.42) 


where [x] q := ( q x — q~ x )/(q — q _1 ). We shall here consider the highest weight 
representation with s + |0) = 0, |0) = s |0) carried by Ti fl = V s , the so-called Verma 

module, corresponding to a complex spin representation parameter s G C: 


s„ = 


sj = 


s„ = 


J2( s ~ k )\ k ) ( k u 


k =0 


E 

k =0 

00 

E 

k =0 


sin (k + 1)77 


sm rj 

sin(2s — k)rj 
sin 77 


\ k ) (A; + 1| 


(2.43) 


\k + l)(k\. 


We stress that the operator only exists within a representation, and not as an 
element of the U q {s\ 2 ) where only q ±sZ exist. Note that only for half-integer spin 
s G ±Z+ = {0, |, 1, |..}, the module V s is truncated, for any 77 , to a common finite- 
( 2 s + l)-dimensional irrep, since | 2 s) then becomes the lowest weight state, s” | 2 s) = 0 . 
The module V s is also truncated to a hnite-m-dimensional one, for any s, when q is a 
generic m-th root of unity, i.e., 77 = 7 r l/m,l,m G Z + . For generic values of parameters 


77 , s, the module V s carries an infinite-dimensional irrep and Lax operator (2.40) is non- 
Hermitian, as in general (s“)f 7 ^ s+, (s^)f 7 ^ s^. 

The matrix R(<p>) satishes the following useful relations: it is symmetric under 
transposition, its derivative at (p = 0 is proportional to hamiltonian interaction, and it 
has a simple inverse proportional to R(—ip): 

RW t = R{v), (2.44) 

d^R{tp) |^=o= |(/i + cosj?l), (2.45) 

R(<p)R(—<p) = (sin 2 77 — sin 2 </?)l. (2.46) 
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Thus, expanding the RLL relation (2.38) for ipip = ip ±8 to first order in 5 one obtains 
a differential form of YBE, or the so-called Sutherland relation 


where 


[hi i2 ,L 1 ((p,s)L 2 ((p,s)] = —2 sin g Li(<p, s)L 2 (<p, s) - Li(<p, s)L 2 (</?, s) (2.47) 


L {ip, s ) = <9 v L(<p, s) = cos cos{i]s z s ) <g) cr° — sin p sin(^s^) ® a z . (2.48) 


This relation (|2.47|) is sometimes also referred to as local operator divergence condition. 

(2.49) 


Let us denote by 

L T (<A s) = s ) ® ( a ' 


a\T 


a£j 


the (partial) transposition with respect to the physical space, noting L T {p, s) = L {<p, s ). 
Sutherland relation transforms under partial transposition in physical spaces to 

[hi, 2 , Li (^,s)L^(<p,s)] = 2sin77 (%(</?, s)L^{p, s) - Lf (<p, s)Z 2 {p, s)) . (2.50) 

We can think of L t (7t — ip, s) also as the Lax matrix corresponding to the transposed, 
lowest-weight representation Vj of U g (s [ 2 ) which has exactly the canonical form (2.40) 
with spin operators transformed under the following algebra-(2.42 (-preserving canonical 
transformation 

„± 


s s —» 


s 7 S s > 


(2.51) 


which is just the spin reversal in auxiliary space. The Lax operator (2.40) is invariant 
under the combined spin reversal in the physical and the auxiliary space 

cr x L(<p, s)a x = L r (7r — ip, s). (2.52) 

ft will turn useful to study also the product of complex spin representations Vj ® Vt, 
for some s,t G C. Namely, defining a double Lax matrix as the following operator over 
the tensor product End(77 a ®77b®77p) with a pair of auxiliary spaces carrying irreducible 
representations of U q (s 1 2 ) with spin parameters s,t G C, "H a = Vj, Rh = V*, and the 
corresponding spectral parameters <p, id G C, 

L x {ip, id, s, t ) = L l x (ip, s)L b>x (#, t), (2.53) 

L x (ip, id, s, t ) = d 5 (L ^ x (ip + 5, s)L b , x (tf - S, t)) s=Q 

= s ) L b t) ~ L)( x (<p, s)d § L htX {id, t ), (2.54) 

where L aiX = L Q ®l b ®o-“, Lf x = L"®l b ®(<x“) T , L b)X = E ae y 

we find again the corresponding Sutherland relation 

[hi 2 , IL 1 IL 2 ] = 2 sin g (IL 1 IL 2 — 1 L^ 1 L 2 


(2.55) 

This identity can be proven directly using the Leibniz rule and Sutherland relations 
(2.47 2.50), or it can be again derived by differentiating YBE for the triple |J(j Vi <g) Vi <S> 

(V®v,), 

Rip (^1 — <5 2 )lLi(<p + ^i; id — S\, s, t)h 2 (ip + 82 , id — 82 , s , t) (2.56) 

= Li(<p + 82 , id — 82 , s, t) 1 ^ 2 ( 9 ? + $ 1 , id — 81 , s, t)Ri, 2 ( 8 i — 82 )- (2.57) 

|| Note that the physical spin space carries the fundamental representation of C/ g (sI 2 ), 7f p = Vi. 
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For a notational convenience, we are using double-strike-roman fonts to designate 
operators which act non-trivially over a tensor product of a pair of auxiliary spaces 
7-La, < 8 > Tib- 


2 . 4 . Matrix product solution - Lax operator method 

Sutherland relations can be straightforwardly facilitated to solve/satisfy defining 
relations ( 2.9[2.10 ) for the amplitude operator. This idea has first been implemented in 
Ref. [3U], even though the corresponding Lax structure and Yang-Baxter equation have 
been identified only later in Ref. m 


Writing (2.50) for a pair of physical sites (x,x + 1), multiplying with Lf ■ • • 

L) from the right, and summing over x — 1 ,..., n, one 


from the left and with Lj +1 


obtains a telescopic series yielding 

[H, L[L^-.-L^] = 2sinr / (L 1 L^--L^ 
Making an ansatz 


Lf • • ' Ln-iL, 


- I n 


sin n (<^ + 7]s) 


( 0 |L^-..L^| 0 ) 


(2.58) 

(2.59) 


one sees that since 

(0| L = (a 0 cos p cos ps — a z sin p sin 77 s) (0| , 
L |0) = |0) (a 0 cos p cos 7js — a z sin ip sin 77 s), 


Q n satisfies (2.9) if 


7 r i£ 

p=~, tan r]s = — -. 

2 2 sm rj 


(2.60) 


(2.61) 


The second condition (2.10) of Lemma 1 is satisfied as well due to normalisation of 


the ansatz (2.59) and the lowest weight nature of representation. Equivalently, since 
[h, a z ® a z ] = 0 one can use another gauge and take a twisted Lax operator L T cr z again 
solving the Sutherland equation. Therefore, another ansatz 

Vn = -y <0| LfLf 

sin ra (^ + ?7s) x 1 1 2 


■i£|o> (^ Z ) C 


(2.62) 


solves the same defining equations (2.9|2.10), and provides the same NESS density 


operator (2.11) according to Lemma 1, if 

ic 

P = 0, cot 7]S = -. 

2 sm 77 


(2.63) 


Of course, the ansatz (|2.59|) provides just an alternative formulation of MPA (2.22) with 

cos(t 7S z ) 


the matrices identified as 


7 T 


A 0 = (sec rjs)L ( —, s) = 


cos 77 s 


A ± = (sec„ S )L = 


tan 77 s = 


1 £ 


2 sin 77 


(2.64) 
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Or, alternatively, picking representation (2.62|2.63): 


A 0 = (cscr/s)L°(0, s) = 


sin(rySg) 
sin r)s 


A± = (cscr]s)L T (0, s) = . \ sf, 

sin (rjs) 


cot rjs = 


1 £ 


2 sin rj 


(2.65) 


Both cases reproduce exactly the result (2.36), up to a gauge transformation | k) —> 
Ck \k ), (k\ —* cf l (k | which does not affect the MPA. For instance, Ck can even be chosen 
to make all transition amplitudes of A^, A linear in dissipation e as in Ref. pi. 

We note that, as a consequence of the spin-reversal symmetry of the Lax operator 


(2.52), the non-transposed, highest-weight Lax operator at (2.61) reproduces the lower- 


triangular amplitude operator kl' n (2.20) for the reverse driving (2.21) 

1 


O' = 


(0| LiLg 


L n |0> . 


sin n (<y3 + gs ) 

2.5. Matrix product solution - the case of isotropic (XXX) chain 


( 2 . 66 ) 


Let us here briefly list the result (appearing e.g. in Ref. |77j ) for SU{ 2) symmetric 
XXX chain with A = 1 which correspond to the limit 7 ] —> 0, properly normalised 
when needed, of the results derived in the previous subsections. The spectral parameter 
we set now as A = (p/p, so the A-matrix and the Lax operator read, respectively, 


A(A) <- 

— lim -R(<p) = 

r/->0 T] 

L(A, s) <- 

— lim —L(o?, s) 

r?->0 Tj 

l(s++s7), 

-|( s ^- s 7) ; s s) ; 


(2.67) 


Ala + 


S+ Al a - 


— A1 + S s • cr, 


over 'H p 0'Hp and s^, have become standard generators of angular-momentum algebra 

3(2 for a complex spin s 


= 


sj = 


S s = 


x>-*)i*><*i« 

k =0 
00 

+ 1 ) ifc)(/c+1|, 


( 2 . 68 ) 


fc =0 

OO 


^(2s-fc)|fc + l)(fc|, 


k =0 


which is genuinely infinite-dimensional, unless s G \7L + . The MPA (2.22) for NESS 


amplitude operator has now generating matrices, obtained from the case (2.65[), as 


qZ qlL 

An =X A+ = 


2i 

s = —. 

£ 


(2.69) 
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2.6. Generalizations of boundary driving 

The solution of NESS for XXZ model presented so far refers to an extremely 
minimalistic boundary condition with a pure source and pure sink of equal rates at 
each end. Using generalisations of our approach we can expand the integrable boundary 
conditions in three different directions: (i) allowing arbitrary left-right asymmetry in the 
source-sink rates combined with additional arbitrary magnetic fields at the boundary 
sites, (ii) allowing the source and sink jump operators to act with respect to non¬ 
parallel axes (like z-axis in our previous example), so the two target states cannot be 
chosen mutually diagonal, (iii) allowing four different rates at the two boundaries for two 
independent in&out processes at each end, but only perturbatively in the system-bath 
coupling constant. We shall describe these developments in some detail in the three 
paragraphs below. 


2.6.1. Left-right asymmetry and combined coherent-incoherent driving. Here we are 
considering the quantum master equation with a combination of asymmetric coherent 
and incoherent driving [59j. The former is provided by adding an arbitrary magnetic 
field at the chain ends 

n—1 

H h = ^2 h x,x+i + b L al + b R a z = H + b L al + b R a z n , (2.70) 


X=1 


while for the latter we allow arbitrary rates of the source and the sink 

L\ = y/lfaf, Lo = v^Tr^, 

so the total Liovillian generator reads C = —iadi/ b + Y L T> + + Y R V -. 


(2.71) 


We start by making the following ansatz for NESS 

n 

Roc = (o, 0| ll L x {<p, 6, s, t ) |0, 0) K® n = (0, 0| (L(p, 0, s, t)K)® n |0, 0) (2.72) 


X=1 


where L x is the double Lax operator (2.53) over the tensor product of a pair of Verrna 
modules "H a <£> Hb = Vj <8) V t and 

T/2 o 


K = K( X ) : = 


X 


MG 


0 X 

is a magnetization-shift matrix satisfying, for any y G M + , 
[h, K ® K] = 0. 


(2.73) 


(2.74) 


Using Sutherland relation (|2.55|) and Eq. (]2.74|) one again implements the telescopic 
series 


to find for the commutator with the entire Hamiltonian ( | 2 . 70 [ ): 
[H h , (L/t)® n ] = (2(smr))lK + b L [a z ,hI<]) ® (LA') 

- (LA') lg>n_1 0 ^2(sin rjfhK - b R [a z ,LK] 
Hence the sufficient condition for the steady state Lindblad equation 

i[/7b,-Roo] = Y L V a +(Roc) + YrD- (Roc) 


-l 


(2.75) 


(2.76) 







Matrix product solutions of boundary driven quantum chains 18 

to hold is to satisfy a pair of boundary equations on A a , <E> Ab <E> 7-L p \ 

(0,01 2i(sinr/)L K + Yifb a +{fLK) — i& L [cx z ,LA']) = 0, 

(2i(sin 7 /)L/i + r R P ff - (LA') - i& R [a z , LA']) |0, 0) = 0. (2.77) 

For convenience of calculations one may write the local left and right dissipators as 
explicit matrix maps over End('Ap) 


V n 


v n 


—2a —b 
—c 2 a 


(2.78) 


a b\ f 2d —b 
c d J y —c —2d 

Let us write the two representation parameters as s,t € C and the corresponding 
generators of Ugisi-f) as 


s ± =af 




s = s: 


Ik. t± = 1 


sf, t = l a 


L = 


(2.80) 


s z , (2.79) 

hence 

sm((p—r)s) sin(i?— rjt) + (sin 2 ?/)s + t + (sm(ip—rjs)t~ + sin(i?+ 77 t)s + ) sin r) 

(sin( , d —?/t)s _ + sm((p+rjs)t + ) sin rj sm(ip+rjs) sin($+?/t) + (sin 2 ? 7 )s~t^ 

~ / sin(</? — 77s — d + 77t) (cos( , d+7/t)s + — cos(<^ — ??s)t~) sin 77 

y (cos(t 9 — ?/t)s _ — cos(y>+?/s)t + ) sin 77 sin(</? + 77s — $ — 771) 

Noting the identities 

s | 0 , 0 ) = s | 0 , 0 ), s + | 0 , 0 ) = 0 , s” | 0 , 0 ) = (sin( 2 ? 7 s)/sin 77 ) 11 , 0 ), 

( 0 , 01 s = s ( 0 , 0 | ( 0 , 0 | s + = ( 1 , 0 |, ( 0 , 01 s~ = 0 , 

t|0,0) =t |0,0), t + |0,0) = 0, t |0,0) = (sin( 27 7 t)/sin 77 ) |0,1), 

< 0 , 0 | t = t( 0 , 0 |, ( 0 , 0 |t+= ( 0 , 1 |, ( 0 , 0 |t" = 0 , (2.81) 


the boundary equations (2.77) amount to two sets of 2 x 2 equations (components 
in End(Ap)), where only 5 are independent for 5 unknown parameters s,t,ip, t},x, 
specifically: 


tan(<£> — 77 s) = 
tan(y> + 77 s) = 
tan(d — 77 1) = 
tan($ + rjt) = 


2 i sin 77 
T l - 2 i 6 L ’ 

2 i sin 77 
Tr + 2 i£> R 7 
2 i sin 77 

Tl + 2i&L ’ 
2 i sin 77 


X 


r R — 2 i 6 R ' 

2 sin(i? + rjt — ip — rjs) 


r L sin(< 7 ? — ijs) sin(-d — rjt) 


(2.82) 

(2.83) 

(2.84) 

(2.85) 

( 2 . 86 ) 


sin(d — rjt — ip + rjs) r R sin(^ + rjs) sin(d + rjt) ’ 
while the other 3 equations reduce to identities. One then easily finds an explicit solution 


tp = d= -( Zl - Z R ), 


_ 1 


rjs — rjt — -(z l + 2 R ), 


(2.87) 

( 2 . 88 ) 


X = X = 


(Kl -^-sm 2 rjy + ribi 


R 


R 

Ur2 
^4 i L 


R 

bt 


Sill' 


: <?) 2 + r Ibl' 


(2.89) 
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where 


1 2i sin rj 

zl ■= t arctan --—. 

1 r L - 2i b L 


1 2 i sin r] 

Zn - arctan —-—. 

i r R + 2 i 6 R 


Note that such NESS is again of Cholesky form, namely defining 


^n(<AS,X) := (°l L l(^, S ) L f(^, S )--- L »(^, S )l°) 


X 


1/4 

o x 


0 

-1/4 


(2.90) 


, (2.91) 


which is a lower-triangular matrix, one can write the non-normalized density operator 
of NESS, since [£)(<£>, s, x)]t = [fl(<^, s, x)] T , as 

Roc = n n (^,s,x)[n„(v7,s,x)] t - (2-92) 


In the isotropic case A = 1 the NESS solution, after writing A = <p/r) and taking 
the limit 77 —>• 0 , can be written compactly as 

Roc A2n(A, s, x) [fl n (A, s, x)]^, 

0,(A,s,x) = {0|L r (A,s)®“|0> A-(^)®”, (2.93) 


with the spectral, representation and magnetisation parameters, respectively 


r L - 2 i 6 L r R + 2 i 6 R ’ 

i i 

r L - 2 i&L r R + 2 i 6 R ’ 
gj + 4 & r) F l 

( r L + 4^)r R ' 


(2.94) 


2.6.2. SU (2) — twisted boundary driving. In the isotropic case A = 1 we shall now 
further exploit the 577(2) symmetry of the problem to map the solution ( 2.93|2.94 ) to 
the NESS for a more general class of dissipators (essentially following Refs. [4DJ [55]). 
We start by elaborating on rotation symmetry of the double Lax operator L entering 
the boundary conditions (2.77). We construct a pseudo-representation of the rotation 
group over "H a ® ®R P = Vj <E> V t 


u (axis of rotation) we define 


Vi, namely choosing an angle 9 and a unit vector 


U(6 I , u ) = exp ^i9u ■ ® 1 2 + t <g) t 2 + l a <E> lb ® — J J (2.95) 

= U s (0, u) ® U t (0, u) ® U(9, u), (2.96) 

where U s (0, u) = exp(idn • s s ), U(9,u) = exp(ifw ■ a), formally implementing j^jj the full 
SL( 2) symmetry of the non-Hermitian double Lax operator 

U(0, u)LU(-0, u) = L, U(0, u)LU(-0, u) = L. (2.97) 


One perhaps needs to stress that the operator U s may not exist as an element of End(R a ) but its 
action on the highest weight state U s |0), (0| U s is well-defined and computable, and that is all what 
we need here. 
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The boundary equations (2.77), substituting RHSs of (2.97) by the corresponding LHSs, 


and using factorisation (2.96), and requiring \ = 1 so that [U(0, u), K] = 0, map to 

«0| U s <g) (0| U t ) (-2i(sin77)L + T L V UrT+uf ( L) - ib L [Ua z U\ L]) = 0, 

'2i(sin7?)L + r R V Ua - ut (h) - ib R [Ua z U\h]) (Uj |0) ® Xj\ |0» = 0. (2.98) 

Note that the two formal SL( 2) transformations for two, left and right boundary 
conditions can be independent, say U(6 , l,ml) and U(0 R , w R ), and without loss of 
generality we may chose the axes of rotation to lie in the x — y plane, ur/r = 


(sin 0 l/r, — cos 0 l/r, 0). Thus the Eqs. (2.98), in combination with 577(2) invariant 


bulk condition (2.75) [for K = 1 2 ] and noting that t = s, implies that the density 
matrix 


fi 0O = a“ i “(A,»,x)K w “U,»,x)] t , 
W =(0|U.(9 l,5l), 

IV’r) = U s (—0 r,ur) |0), 


(2.99) 


( 2 . 100 ) 


with A and s determined from the first two lines of (2.94), represents an exact NESS of 


the Lindbladian dynamics with the twisted jump operators 
L twisted = U(0 l ,u l )<t+U(-0 l ,u l ) 


( 2 . 101 ) 


3 - i^l 


- (cos Or cos 0 l — i sin 0 L , cos Or sin 0 L + i cos 0 l, — sin 9r) ■ 3 \, 


- twisted 


oW R 


— U(0R,UR)a n U(—O r, Ur) — 

(cos Or cos 0r + i sin 0 R , cos Or sin 0 R - i cos 0 R , - sin Or) ■ a n , 


and for the Hamiltonian with twisted boundary fields: 


n—1 


= ^2 h XiX+1 + 6 L (sin Or cos 0 l , sin Or sin 0 L , cos 0 L ) ■ a x 


X=1 


+ 6 R (sin Or cos 0r, sin Or sin 0 R , cos Or) ■ a n . (2.102) 


The states (0 R | and |0 R ), (2.100), are just the SU( 2) coherent states over the complex 
spin Verma module, namely 


(0 L | = (0| exp(0 L s+ - ip L s s ) = [cos 0 r\ 2s ^ ^ 


k=0 
oo 


-tany) e- ifc0L 


0 


R 




\ k ) , 


|0r) = exp(0 R s+ - 0 R S S ) |0) = I cos 6 > r | 2s ^ ( - tan 

k=0 

where the complex coherent-state parameters read 
0L = -ye-^, 0 R =ye- i( K 
The expansions are derived using a well known SU(2) disentangling formula: 
exp(0e i *s + - 0e~ irt> s") = exp(-s”e“ i0 tan 0) exp(2s z log [cos 0|) exp(s+e^ tan 0). (2.105) 


(2.103) 


(2.104) 
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Matrix product solutions of boundary driven quantum chains 
Note that solvability condition (y = 1) in this case 

(Fl + 4 /jDFr = (rj + 4 / 4 )r L . (2.106) 

still allows some left-right asymmetry, in which case our solution goes beyond what has 
been discussed in Ref. [4Q S 63]. Due to rotational symmetry, we could of course without 
loss of generality fix three out of four angles 0 l = 0l = $r = 0 and keep only the relative 
twisting angle 0 r between the source and sink measurement axes [63]. Nevertheless, it is 
perhaps illuminating to write out the full rotationally symmetric parametrisation of the 
twisted boundary driven solution for possible further generalisations or deformations. 

As we have seen, general analytic solutions for twisted boundary driving are limited 
to the isotropic case A = 1 only. However, interesting nontrivial properties of the spin 
current under twisted driving have been observed in the anisotropic case A / 1 in 
Ref. [62] by exact analysis of short chains. The question of non-equilibrium integrability 
in such a case remains open. 


2.6.3. Perturbative driving with source&sink processes on each end. The explicit forms 
of NESS of boundary driven XXZ chains that we discussed so far were all characterized 
with simple, ultra local, rank-one dissipators, which can be considered as a source and 
a sink of spin excitations with respect to some measurement bases. I 11 the classical 
integrable locally interacting Markov chains, e.g. ASEP with open boundaries [8], one 
can analytically solve more complicated boundary conditions with in and out incoherent 
processes on each side. I 11 the framework of Lindblad equation, these would be described 
by four Lindblad channels, with four non-negative rates T^ R > 0: 


+ ^-+ 

1 ? 




Li = 

generating the Liouvilllian C = 


Lo — 




1 5 


— 


' r K- 


u = 


r R a-, (2.107) 


-iad//b + 'Yhu=i'i- ) L il in terms of Hamiltonian (2.70) 


and canonical dissipators (2.6). In analogy to ASEP, one may hope that the NESS 
can be written with an ansatz generalizing (2.72), R^ = (Tl| (LJi )® n |Tr), where 
|T l /r) G 'H a ®PLb are some free auxiliary states. However, a straightforward calculation 
performed by the author showed that the above ansatz is insufficient, i.e. there is 
generally no solution for |Tl/r) and parameters <p, 6, s, t , x■ Therefore, finding an exact 
solution of NESS for XXZ chain driven by four channel boundary dissipation remains 
a challenging open problem. 

What one can do instead is to solve the problem perturbatively if all the coupling 
and driving rates are uniformly small (see Ref. [H]). Writing 

b-L/K = £Rl/r, (2.108) 


p± 

1 L/R 




where £ is a formal small parameter and expressing an un-normalized density operator 
as a power series 


R x = ^(i£)V P) 

p =0 


(2.109) 





Matrix product solutions of boundary driven quantum chains 


22 


we get a recurrence relation connecting subsequent orders (where H below now denotes 


the free-boundary XXZ Hamiltonian (2.2)) 

(ad H> (0) =0, 

(ad H)p^ =—V(p < ' p ~ 1 ' 1 ), p= 1,2..., where 
V = 7l ^ 0 -+ + 7l Arp + 7 r Ar+ + 7 R^a- - i/JR ad a 7 n - i// L ad < 

From the uniqueness of NESS it follows that each term p ^ in the formal series expansion 


( 2 . 110 ) 

( 2 . 111 ) 


(2.109) should also be determined uniquely. This means that even though at each fixed 
order p, the solution of Eq. (2.111) p^ is undetermined up to addition of an arbitrary 


element from the kernel of ad H (operator which commutes with H ), there is always 
a unique element p H such that 'D(p^) is in the image of ad H (it is Hilbert-Schmidt 
orthogonal tr{X^I?(p^)} = 0 to all operators X which commute with H ), so that the 
equation in the next order p + 1 can have a solution. Note that the map ad H is self- 
adjoint w.r.t. Hilbert-Schmidt inner product hence the orthogonal complement of its 
image is its kernel. 

In the leading nontrivial order, the perturbative solution of NESS can be encoded 
compactly in terms of the Z-operator [69], a strictly upper triangular matrix Z G 
End (fH® n ) which satisfies a remarkable conservation law property 

[H,Z} = -a\ + al. (2.112) 

The operator whose time-derivative is composed of local operators at the chain 
boundaries has been termed as almost conserved [35] and provides a useful tool to study 
the thermodynamic limit of quantum transport via the Lieb-Robinson bounds [50]. Our 


operator Z can be expressed in terms of a derivative of amplitude operator (2.22) w.r.t. 


noise strength, or highest-weight transfer matrix w.r.t. representation parameter at 
s = 0 [70] 

1 


z '^d £ k} n \ £= Q 


cl 


7T 

2 ,S 


|o> U 


(2.113) 


2r/ sin r] 

and has a simple explicit MPA representation [69], related to ( 2.31[2.36 ) at e = 0 in an 
extended auxiliary space with a split-vacuum state T-L' & = lsp{|L), |R), |1), |2)...} 

z= E (2.114) 

ai ,0} 


Ao = L)(L + |R)(R| + E cos (H \k)(k \, 

k= 1 
oo 

A+ = |L)(1| + sin(/cr]) | k) (k + 11 , 

k=1 

OO 

A'_ = 11 )(R| - ^2 sin ((l + 1 )v) \k + l)(l| ■ 


k =1 


To first order, up to 0(e 2 ), the following simple ansatz 

p(0) = ir(x )®n j p(l) = C (Z - Zt) K(x) m 


(2.115) 
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solves Eqs. ( |2.111 ), which results in a single condition 
C[H,Z-Z t] = 


-ltK?V 0 


K^-K^V-Kr 


- TS~ 1 


-^KpV.K^K, 


„ V„-K„ (2.116) 


imposing a vanishing linear combination of 1, <jj, o 7 n . Requiring the coefficients to vanish, 
results in a system of equations for y, £ with the unique solution: 

1 7l + 7r + 7l + 7r 




7r 


C=X( 


-=\(7l7r -7l7r)- 


(2.117) 


7l + 7r ' J 2 (7l + 7r)(7l + 7r) 

It is worth to remark that these leading order terms of NESS do not depend on coherent 
driving parameters /xl/r- Note that (y — 1)/(y +1) gives the net magnetization tr(<j^p 00 ) 
in NESS to leading order, while e( is essentially the spin current tr[(icr+cr7 + i + h.c.)poo] 
within the first order. In order to obtain the spatial modulation of the magnetization 
density profile, one needs to compute the second order p = 2. With the tools at hand, 
this is only explicitly possible in the case of symmetric incoherent driving 

7 f = i(l± M ), Tr = |(1 =F A), «./r = 0, (2.118) 

where the solution reads (for a simple proof see Ref. [69] ) 

,2 


p(°) = l. p W = p(Z-Z f ), p 


( 2 ) _ f£ 


= y(2-Z f ) 2 -! (2.119) 


In order to obtain a closed form expression for p^ for general driving parameters one 
probably needs to include second derivative of highest-weight transfer matrix w.r.t. s 


at s = 0 (extending (2.113)) However, this has not been explicitly demonstrated yet. 

One can also study asymptotics for a large coupling parameter e —> oo, in the so- 
called quantum Zeno regime, by writing a formal operator valued expansion of NESS in 
1/s, Poo = Yl^Lo £ ~ P ■ ^ ie case ftdly anisotropic Heisenberg spin 1/2 chain 
(.XYZ model) a remarkable effect has been demonstrated [(55], namely engineering 
a transitions from equilibrium-like (uncorrelated) to genuine nonequilibrium (strongly 
correlated) steady state by applying local magnetic fields to spins near the boundary (at 
sites x = 2 and x — n — 1). It is possible also to derive explicit asymptotic expansions 
in some other (large) parameters of the model, say in (spatially) modulated external 
magnetic field or anisotropy A. Even in the generic, non-integrable situation of XXZ 
spin 1/2 chain with spatially modulated interactions an explicit asymptotic expression 
for the spin current has been derived [IB], exhibiting a strong rectification effect upon 
switching the direction of coherent driving in the presence of incoherent driving (see 
also Ref. [SjTJ). 


3. Nonequilibrium partition function and computation of observables 

Here we shall elaborate on computation of physical observables using the standard 
‘transfer-matrix’ technique. For concrete examples, we work out spin-density profiles, 
spin currents, and two point spin-spin correlations. For most of our discussion we allow 


the NESS to be of the most general non-perturbative form (2.72), or equivalently (2.92), 


with parameters (2.87 2.89), as driven by a general combination of asymmetric coherent 
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and incoherent boundary couplings ( 2.70||2.71 ). Only for some very specific explicit 
calculations for the XXX chain discussed in the second part of this section we shall 
assume purely incoherent and left-right symmetric driving. In parts this section follows 
Refs. [701 63] . while some essential results are new and presented here for the first time. 


3.1. The case of generic XXZ chain 

Let us define the nonequilibrium partition function as the trace of the un-normalized 


density operator of NESS (expressed as (2.72)) 


Z n = tvR oc = (0,0\T n 10,0), 

where T = T(<^, s, x) € End(7f a <8 TL\f) is the transfer-matrix 
T = tr p (L/i) = 2L n cosh k + 2L Z sinli k, 
writing the magnetisation parameter x in terms of another parameter k as 


(3.1) 

(3.2) 

(3.3) 

Any physical observable A e End("Hp n ) can be expressed in terms of a linear 
combination A = s f2,a a oPQ L of tensor products = a ai <E> cr" 2 ® ■•■cr" n , a x e J. 
For each product operator 0& its NESS expectation value can again be conveniently 
expressed in terms of a matrix product 


X = e 


Jl K, 


(^a) tr(OaPoo) 


tr {OJl 0 


= Z~ v (0, 0| Y ai V Q2 • ■ ■ V Q " |0,0), (3.4) 


trRoo 

where V" = tr p (L/icr Q ), or, explicitly: 

V ± = e ±K lY, V° = T, V z = 2L Z cosh k + 2L° sinh k. (3.5) 

The operators V" shall sometimes be referred to as vertex operators. Since the generators 


s, s ± , t, expressing the physical components of the double Lax operator (]2.80|) are 


tridiagonal matrices, one needs to consider, for a chain of n sites, only auxiliary basis 


states |A;) up to k < n/ 2, and consequently the above expression (3.4) can be evaluated 
efficiently within 0(n 3 ) computer operations even without any further insight. We will 
anyway show bellow that expectation values of many observables in many situations 
can be evaluated fully analytically. 

We shall particularly focus on three kinds of observables in NESS. The simplest 
and perhaps the most important one is the spin current 


jx,x +1 — i(°"x" Cr a-+l TrTr+l) 


(3.6) 


which satisfies the continuity equation, or local conservation law for spin density <r z 


—al = i = -4j„ +1 


(3.7) 


We shall use (3.6) as the spin-current operator in the following, although we note a 


trivial factor of 4 which needs to be taken into account when comparing to physical 
units. The above identity also implies that the expectation value of the current should 
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be site independent in the steady state (j x , x + 1 ) = {j 1 , 2 ) = : J■ Similarly as in the 
case of ASEP, the current in NESS can be computed solely in terms of nonequilibrium 
partition function Z n . To see that, we express the asymmetry parameter y from the 
defining equations (2.82 2 . 86 ) as a simple function of four free complex variables <p, it, s, t 

= + (3.8) 

V sm(<p — it + ps — rjt) 

We note that arbitrary branch of the square-root can be chosen as it only affects the 
sign of the transfer matrix T, whose explicit form we read from (2.80): 

T = (sin 2 ? 7 )x 1 / 2 s + t + + (sin 2 7 /)y _ 1 // 2 s~t _ 

+ x 1 / 2 sin(<p — 77 s) sin(i? — 77 1 ) + y^ 1 / 2 sin(<p + 77 s) sin($ + 77 1 ), (3.9) 

or the sign of the un-normalized density operator for odd n, but not the observables 
(3.4) themselves. One notes that a similar expression is obtained for the commutator 
of the off-diagonal elements of double-Lax operator (csc 77 )[L _ , L + ] = (csc 77 )[V + , V - ] = 
(sin 2 77 ) sin(<p — it + 77 s — ? 7 t)s + t + + (sin 2 77 ) sin(d — p + 77 s — 77 t)s _ t _ — sin( 2 ? 7 s) sin(d + 
r/t) sin(d — ? 7 t) + sin( 2 ? 7 t) sin(<p + 77 s) sin(<p — 77 s). 

Furthermore, let us identify a particularly important, diagonal subspace of the 
product auxiliary space /C = lsp{|/c, k ), k = 0,1, 2 ...} C "H a < 8 > Rb, where a compact 
Dirac notation | k,l) = | k) ® |Z) is used. We note that any operator valued function of 
/(s — t) on /C evaluates as f(s — t). Henceforth, elementary trigonometry results in the 
following very useful relation for the orthogonal projection on /C 


P([¥+,¥~] -iCT) = ([¥+,¥“ 


= 0, 


OO 


P :=l^\k,k)(k,k\, 

k =0 


where 


(3.10) 


(3.11) 


C(< p, it,s,t ) = (sin / 7 )\/sin(<p — it — 77 s + 771 ) sin(<p — it + 77 s — rjt). 

Now using the definitions ( 3TpH ), together with the facts [T,P] = 0, P10,0) = |0, 0), 
one shortly arrives at a compact expression for the steady-state spin current 
.Z n -\ 


J = c- 


Zn 


(3.12) 


which is similar to the expression of a particle current in the classical simple exclusion 
processes [B] • Note that the parameter ( can be conveniently expressed also in terms of 
physical driving parameters, via (2.82||2.85|), namely 


c = 


(riJn ) 1 / 2 


sin 2 rj 


(3.13) 


((jp ->i- sin 2 o) 2 + r 2 L 6£)‘ /4 ((hi - b i- si “h) 2 + r|6 |)' ,,r 

It the simplest case of symmetric incoherent driving, Tl = Tr = e, = &r = 0, it 
amounts to ( = esin 2 77/||£ 2 — sin 2 77 1. 

Other simple and interesting physical observables that we consider in some detail 
are the spin-density and connected transverse spin-spin correlation function of NESS 

M x = (a z x ) = Z~ l (0, 0| T x_1 ¥ z T" _a: |0, 0), (3.14) 

C x , y = ( <J z x a z v ) - ( a z x )(a z v ) = Z~ 1 <0,0| |0,0) - M x M y , 
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where in the last line we have assumed, without loss of generality, that x < y. Since 
V z also conserves the spin-difference s — t and hence leaves the diagonal subspace /C 
invariant i.e., [V Z ,P] = 0, one can identify the diagonal subspace with the auxiliary 
subspace /C "H a via the mapping \k,k) | k), and write the diagonally projected 

transfer matrices T := T|^, V :=V z |yc as 


T = V 00 + V 11 , V = V 00 - V 11 , (3.15) 

00 

V 00 = x _1 ^2 ( I sin ((^ + l) 7 ?) 1 12 I k)(k + 1| + I sin(</? - 7](s - k))\ 2 \k)(k\^j, 

k =0 

OO 

V 11 = x£(l sin((2s — k)rf) f \k + l)(fc| + | sin(<^ + p{s — k ))| 2 \k)(k\^j. 
k =o 


We explicitly used the complex conjugation property rjt = ffs, $ = (p, which one has 
for physical values of the driving parameters (2.87). In terms of the projected transfer 
matrices, the nonequilibrium partition function and the transverse spin observables read 


= (0|T"|0>, 

M x = Z- 1 (0| 

T :r_1 VT y ' 


C —7 
x,y ~ z -’n 


(3.16) 

|0), (3.17) 

VT" -1 ' |0) - M X M V . (3.18) 

For the massless XXZ model, |A| < 1, one can always approximate 77 = arccos A e 
R to an arbitrary accuracy, for fixed n, with a rational q = nl/m, with l, m G Z, m > 0. 
This corresponds to q — e lv being a (generically non-primitive) 2m—th root of unity. 
In such a case, the transfer matrix T can be truncated to an m-dimensional sub-space, 
= lsp{|/c), k — 0,1 ..., m— 1}, since the transition between states | m) and |m — 1) is 
forbidden (see the first summation term of V 00 in (3.15)). Hence the transfer and vertex 
matrices T. V, can be replaced, respectively, byroxm matrices, T' = T| %/, V' = V| 


I H' 


Specifically, T x |0) = (T 7 )® |0), \/x, so in TL, n —> 00, observables are essentially given 
by eigenvalue decomposition of T' = U diag(r 1; r 2 ,... r m )U -1 , with eigenvalues ordered 


as |ri| > IT21 > 




T m I, so the steady-state current is ballistic (i.e., n-independent) 

(3.19) 


r, 


Similarly, using the fact that in the eigenbasis of T', the transformed vertex operator 
U _1 V / U have vanishing diagonal elements 0- Wj\ V' 1 4>j) = 0 for T I fjj) = Tj 
(fjl | T 7 = Tj (ipj |, eigenvalue decomposition gives thermodynamically vanishing spin 
density with exponentially damped profiles near the ends (Fig. [3^i), namely 




3 =2 



X — l 


Cd \ I — 


with Cj = (0| ifj) (ifj | V' |V>i) <^i|0). (3.20) 


+ This follows from a simple observation that there exist a gauge transformation | k) —X </>£ | k), 
(7j — > (fk)- 1 ( k\ for appropriate weights <pk 0 such that T' becomes a symmetric and V' an anti¬ 
symmetric matrix, T ,t = T',V' t = —V'. 
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Similarly, it can be shown that far away from the edges, 1 <C x, y <C n, spin-spin 
correlations decay exponentially 

\x-y\ 


C « 

^x,y ~ 


3 =2 


with d = (OlV’i) {'if il V' | i/)j) {tf/j | V' l^i) (V’ilO). (3.21) 


Note that such exponential decay of correlations in the steady state in this regime is 
qualitatively reminiscent of equilibrium behaviour at a finite temperature T > 0, namely 
long-range order is absent, even though p^ is highly non-thermal (non-Gibbsian). 
Let us now work out an explicit example for A = 1/2 = cos(7t/3), m = 3, and 
for symmetric incoherent driving IA = Tr = e, b^ = 6r = 0, implying ip = 0 
and tan qs = 2i sin p/e = i\/3/e. Up to trivial similarity (gauge) transformation 
| k) —» (||e 2 — 3|) 1 | k), {k | —> (|/ 2 — 3|) {k |, the transfer and vertex matrices read 

l^ 2 o . 

+ ^(l + ^ 2 )(9 + e 2 ) I , (3.22) 


r = 2C 




±E 2 

4 C 

0 

-1 


K^ 2 ) 

0 


^(l + e 2 )(9 + e 2 ) 


(3.23) 


The eigenvalues of T' are ri )3 = ^(7 + 3e 2 ± a/81 + 74e 2 + 9e 4 ), r 2 = |(5 + £ 2 ), yielding 
the spin current J = 8^/(7 + 3 e 2 + a/81 + 74 e 2 + 9£ 4 ) (see Fig. [3|d). Spin-prohles (see 
Fig. |3^,) and spin-spin correlation are described by Eqs. ( 3.20|3.21 ) with explicit, but 
lengthy expressions for the coefficients 02,3(5), d 2 3 (e). 

On the other hand, for the massive XXZ model, |A| > 1, parameter r/ is complex, 
namely 7] = h/, with rf = arcosh |A| for A > 1 and rf = arcosh | A| + in for A < — 1. In 
such a case, the tridiagonal transfer matrix T (3.15) is of genuinely infinite dimension, 
with exponentially growing (in state index k ) transition amplitudes. In WGS picture 
[see Eq. (2.37)] the walk that gives a dominating contribution to the partition sum 
Z n , for large even n, is composed of n/2 steps forward (0,1), (1, 2),..., (n/2 — 1, n/2) 
followed by n/2 steps backward (n/2, n/2 — 1),..., (1, 0), 

n/2 

■Z"n — n | sinh((2s — k + 1)t/) sinh(7c?7 / ) | 2 . 


(3.24) 


k= 1 


We note that the contribution from this extremal walk relatively overweights the sum 
of all other contribitions for large n, as it grows faster than any exponential in n, 
so it is super-exponentially larger than contributions of exponentially many (< 3 n ) 
typical terms. The spin-current, for large n, can then be computed by applying (3.12) 
twice, namely {(/J) 2 = Z n jZ n _ 2 = | sinh((2s — n/2 + 1)?/)| 2 | sinh(n7//2)| 2 , yielding 
asymptotically J ~ C > /\ sinh(n7//2)| 2 ~ (\e nv ’\, or (see Fig. for comparison with 
exact numerical results from transfer matrix computation) 

c 


J 


(|A| + a/A 2 - l) r 




(3.25) 
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Similarly, one can compute the spin density Ai x , obtaining a kink profile from the 
dominating walk, namely first nj 2 spins pointing up and last nf 2 spins pointing down, 
M. x = sign((n + l)/2 — x ), where sign(x) = —1, 0,1 for x <, =, > 0, respectively, while 
connected correlations thermodynamically vanish C XjV = 0. In fact, from the dominating 
walk one obtains the leading asymptotics for the entire NESS density operator 

poo - K + ^r) • • • «/2 (T n/ 2 )( CT n/2+l CT n/2 + l) ' ' ' ( a n a n ) 

= |00... Oil... 1) (00... Oil... 1[. (3.26) 

For Enite n, the kink spin density profile attains a hnite width (see e.g., Fig. [3^,), which 
can be quantified to be of order log n [6J. 


3.2. Isotropic (XXX) chain - nonequilibrium partition function. 


At the end, let us turn to perhaps the most interesting case of SU(2) symmetric XXX 
or isotropic Heisenberg chain. Here the double Lax operator have to be defined with 

lim^oo j]^ 2l L(r]X, r)n, s, t ), reading 
(A — s)(fjj — t) + s + t + (A — s)t _ + (/r + t)s H 
(A + s)t+ + (yU — t)s (A + s) (yU + t) + s t 

— lirn^oc 7]~ 2 T(r]\, s ) is again manifestly 
infinite-dimensional, but with amplitudes growing only quadratically with the state 
index k 


respect to a scaling limit L(A ,p,s,t) i — mn, 4oo 


L = 


while the projected transfer operator T(A, s ) <— um ?4O0 


(3.27) 


T = 


OO 

J2(x 1/2 (k + l) 2 1 k + l)(k\ + x~ 1/2 \k - 2s| 2 1 k)(k + 1| + 


k =0 


(3.28) 


+ (X 1/2 \k - s + A| 2 + y 1/2 | k - s - A| 2 ) 

We shall now present a simple scaling argument which allows to analytically compute 
large n asymptotics of the nonequilibrium partition function Z n and consequently the 
spin current, for arbitrary driving parameters, deriving the result announced already 
in Ref. [70] . Let us define a tridiagonal operator T on the space of sequences of 
coefficients ^ = (if 0 , ifi, if 2 ...), with Y^k=o( T if) k \ k ) = T T.T=o V’fc 1^), namely 

C Tip) k = x 1/2 k 2 'ifk-i+X~ 1/2 \k-‘2s\ 2 'if k+1 +(x 1/2 \k-s+X\ 2 +x~ 1/2 \k-s-X\ 2 )'if k . (3.29) 


The partition function (3.16) can be written as Z n = y/f 1 ' 1 where if^ = T n { 1, 0, 0,...). 
We shall however compute large n asymptotics of the entire sequence 'if^. Since T is 
a tridiagonal operator the vector tf^ is supported on exactly n + 1 sites, i.e. if^ = 0 
for all k > n. We thus propose the following scaling a n sat z 

if k %) - F n exp (- nf(k/n )), (3.30) 

where F n is a sequence of real numbers and some smooth (differentiable) function on 
£ G [0,1]. The consistency of the ansatz is demonstrated, and difference and differential 
equations for F n and /(£) are, respectively, derived, from expanding both-sides of local 
scaling relation if ^. n+1 ^ = F n+1 e^ n+1 ^^/( n+1 )) = ( Tif^)k in 1/n, namely 

F n+ ie ( n + 1 )f( k / n )-( k / n )fX k / n ) ~ e n f(X/n) (^ e f {k/n)-K _|_ e -f'(k/n)+K _j_ _|_ g -K\ (3.31) 
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The dependence on parameters A and s can be neglected to leading orders in 1/n, namely 
they yield smaller or comparable correction than neglecting the second derivatives due 
to shifts k —>• k ± 1 or scaling k/n —> k/(n + 1) in the exponentials on the RHS or LHS 
of (3.31), respectively. Introducing a scaling variable £, via k = n£ and dividing by 
F n e n k we finally obtain 


ex p(/(£) - £/'(£)) = 2 £ 2 (cosh(/'(£) - k) + cosh k ). (3.32) 

As RHS does not depend on n, neither must the LHS, i.e., F n+1 /(n 2 F n ) = C , while 
without loss of generality one may fix C = 1 by suitably adjusting /(£) by an additive 
constant. Hence we arrive to 


F n = Fi(n!) 2 , 

and a curiously-looking implicit differential equation for g(£) := /(£) — k£ 
g(€) = 21og^ + log(2cosh^(0 + 2coshft). 

Using a substitution for a new independent variable 

t = -d 


(3.33) 


(3.34) 


(3.35) 


and writing the parametric dependences as gt and £*, the equation (3.34) transforms to 

(]t = 2 log ~ t£t + log(2 cosh t + 2 cosh k). (3.36) 

An explicit differential equation for is obtained by differentiating with respect to t 
and using (3.35) 


2 d & _ , 


sinli t 


(3.37) 


cosh t + cosh k 

This equation can be linearised by substitution y(t) = l/£ t and solved explicitly in 
terms of elliptic integral of the first kind F((j), k ) = dd(l — k 2 sin 2 d) -1 / 2 . Fixing the 
integration constant by incorporating the boundary condition 

= 0, 6^00 = 1, (3.38) 

which expresses the obvious fact that the scaling variable £ has to span the entire interval 
[0,1], one obtains an explicit result 

-l 


6 = 


1 + cosh K 
cosh t + cosh k 


K 


K tanh — + LF —, sech — 


it 




(3.39) 


where K{k) = F{tt/ 2, k) is the complete elliptic integral of the first kind, which together 
with Eq. (3.36) yields the complete scaling profile. From (3.39) we obtain the key 

cosh(«;/2) 


information 


/ (0) = i lim o / t = 21 og if(tanh2(K/2)) , 


(3.40) 


which yields the asymptotics of Z n , and consequently, the spin-current (3.12) 

2 n 


z ~ F e”/<») = F,M) 2 ( co3h(K/2) j “ 

” ” 1 ’’ fif(tanh 2 (K/2 ))) ’ 


f cos1i(k/2) 

= V 2 


k K[ tanh 2 (ft/2)) / n 


c 


y, (3.4i) 















Matrix product solutions of boundary driven quantum chains 


30 


where the AAA-scaled current parameter £ (3.11) reads 

c- 


limn 2 C(??A,7 ]p,s,t) = a/(A~ 

T )—>0 


or in terms of driving parameters 


c = 




p — s + t)(A — p + s — t), (3.42) 


(3.43) 


In the special case of symmetric driving x = 1, k = 0, the scaling profile simplifies 
(noting that K (0) = 7t/2) 


6 = 


2sech(f/2) 


7T 


4 arctantanh(f/4) ’ 
and /(0) = 2 log yielding the spin current [70] 

ttC 


ft = 9t = 2 log ( 2£ t cosh - 


J = 


4n 2 


*6, (3.44) 


(3.45) 


(3.46) 
~ 100 


One finds an excellent agreement of the whole scaling profile 
(k\ T" |0) ~ Fi(n!) 2 exp ( nf(k/n )) 

with numerical iteration of the transfer operator T even for relatively small n (n 
for r l,r, hpR ~ 1), where the undeterminable constant F\ quickly becomes irrelevant 
due to the super-exponential growth (see Fig. [4]) . One can repeat our scaling analysis for 
the transpose of the transfer operator to show the same leading order in 1 /n asymptotics 

(k\ T" |0) ~ (0| T n \k). (3.47) 


Note that the asymptotics of the partition function (3.41) is unique and depends 
only on the asymmetry parameter k and not on the spectral and the representation 
parameters, A ,s separately. In fact, an analogous asymptotics should be obtained, 
following essentially the same derivation, when starting from an arbitrary local state |Z), 
l G {0,1,..., k}, hence giving the scaling of an arbitrary matrix element of T 11 , namely 


(n 


n 2 


(k | T"- ( 10 ~ exp (nf(k/n )), l < k, 


(3.48) 


where F is a constant independent of k,l,n. Using expressions ( 3.46|3.47|3.48 ) one can 
control the asymptotic n —» oo behaviour of any \ T n \^r) if at least one of the 
auxilliary states |'0l/r) G ZC has a finite support , i.e., in can be expanded in finitely 
many basis states \k). 

As we shall see later, the same universal scaling of nonequilibrium partition function 
and the corresponding canonical current Z n _i/Z n oc J oc n~ 2 applies to several other 
models with intrinsic (and undeformed) £77(2) symmetry, such as the nonequilibrium 
Hubbard model or even spin-1 Lai-Sutherland chain. 
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x 



Figure 3. (from [70]) Spin profiles M x = {a x ) at n = 100 (a), and spin currents J 
vs. size n (b), of boundary driven XXZ spin 1/2 chain for A = 3/2 (dashed), A = 1 
(dotted/blue), A = 1/2 (full curves), all for three different incoherent spin source/sink 
rates IY = Tr = e = 1,1/5,1/25 using thick, medium, thin curves, respectively. 
Red full curves show closed-form asymptotic results [see text]: Ai x = cos7r^5yi 
J = 7r 2 £ -1 n -2 for A = 1 in the main panels (a,b), and J oc e - raarcoshA j n (b)-inset. 



Figure 4. The universal scaling profile /(£ = k/n) ~ Mog((fc| T" |0) /(n!) 2 ), 
generating the nonequilibrium partition function, for the boundary driven XXX spin 
1/2 chain with symmetric driving Tl = Tr = 1, 6 l /r = 0 (corresponding to k = 0). 
The points show numerical data for n = 64 (brown points), n = 256 (blue points), and 
n = 1024 (red points), compared to the universal analytical result (3.44) depicted with 
a black curve. 
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In the isotropic case one can also calculate the spin-density profile and spin-spin 
correlations analytically, at least in the case of symmetric (and un-twisted) driving, 
i.e. in the absence of spectral parameters, p = A = 0 and \ = 1. This corresponds to 
driving with Tl = Tr = e, 6 l = — &r = b, yielding the representation (spin) parameter 
2i 

s = -—. (3.49) 


e-2i b 

In such a case, one finds a remarkable algebraic relation between the transfer and vertex 
operators, T = s + t + + + 2st, ¥ z = s + t + — s _ t _ , 

[T, [T, V z ]] + 2{T, V z } = 4(s(s + 1) + fit + l))¥ z , (3.50) 

which can be derived straightforwardly using only SU( 2) commutation relations and 
our complex spin representation. Note that the relation holds even for a tensor product 
of two abstract SU( 2) algebras, where s(s + 1) and fit + 1) have to be replaced by 
the corresponding Casimir operators s~s + + s(s + 1) and t~t + + t(t + 1), respectively. 
Multiplying Eq. (3.50) by (0, 0| TP~ 2 from the left and by T n ^ x ~ l 10, 0) from the right, 
using the definition (3.14) of (adding explicit notation of dependence on the chain 
length n), and noting that t = s, one gets 


= 8Res(s+l)M { fiT 1 2)Zn 2 


_______ , • (3.51) 

'Z-'n 

This is a closed form difference equation for Afi, knowing Z n _\j Z n ~ 7t 2 /(4n 2 ), which 
has a unique solution once we specify the boundary conditions Af 1 )" 1 and . These 
are givien by the following trivially satisfied boundary equations 

(0, 0| (T — V z ) = 2|s| 2 (0, 0|, (T + ¥ z ) |0, 0) = 2|s| 2 |0, 0), (3.52) 

namely multiplying them, respectively, by TF 1 ^ 1 |0,0) and (0, 0| T n_1 , we obtain 


1 + M S n) = 2|s 


,Z, 


n —1 


7 T \S\ 


1 -Mi n) = 21 


, 2 , 


n—1 


z„ 


ft \s\‘ 

2 n 2 


(3.53) 


2 n 2n 2 

We can take TL n —> oo of equations ( 3.51|3.53 ) obtaining the differential equations for 
the scaled spin-density profile 

x — 1 


Mf = 


n 


- A/fW 

- J l X 1 


specihcally 

M"{fi) = -ft 2 M{£)i M( 0 ) = 1 , M(l) = — 1 . 

The bulk and boundary conditions are all correct to order Ofie 2 + 4b 2 )~ 1 n~ 2 ) 
cosine-shaped solution of the spin-density profile 

Mfi) = cos(tt£), 

should be universally valid for any fixed e > 0, in the limit n —y oo, or e e* ~ 


(3.54) 

(3.55) 

. The 

(3.56) 
l/n. 
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We can make a similar computation for the two point spin-spin correlation function 
C Xj y, however here we need to keep the first two leading orders in the 1/n expansion. In 
principle we should now need also 1/n-corrected scaling of the partition function 


Z n ~\ 

Z n 


7 r 


4(n 


a) 


ll + 0(n-*)). 


(3.57) 


It will turn out that the final result for connected correlator C x/y does not depend on the 
value of a as it cancels out from our calculation, so we may leave it as a free, unspecified 
parameter. Nevertheless, numerical simulations suggest clearly that a = 3/4 [70]- We 
start by upgrading the accuracy of the 1-point function A4(£). Expanding Eq. (3.51) 
via (3.54) to 0(n~ 2 ) results in the following differential equation correcting (3.55) 

A4"(£)+7r 2 M(£) = ^ 0 BM(t) + (1 - 2£)AT(£)), 0 := 4(1—ck). (3.58) 

Writing the solution as A4(£) — cos7t£ + n~ 1 A4(f) + 0(n~ 2 ), one finds inhomogeneous 
equation for the first order term 

r 2 


71 


M "(£) + tt -M(£) = — (/3costt£ - tt(1 - 20 sin O , 


(3.59) 


with boundary conditions A4(0) = A4'(0) = A4(l) = A'f(l) = 0, following from further 

(3.60) 


expanding Eq. (3.53), with a unique solution 


7r 


M(g) = - «(1 - £) cos7 t£ + ((1 + /3)£ - 1) sin7 t£) . 


Going next to 2-point function we start with the difference system [following from (3.50) 
for its unconnected part A4x’y = (0,0| , f :c_1 V z T y_x_1 V z T" _1_2/ 10, 0), x < y: 


(n — a 
t(n) 




2 MU + M *h>) + y (^Kf, + = °( 

M?J - MP = C(rr 2 ), - M'p = 0(n- 2 ). 


- 2 \ 


(3.61) 

(3.62) 


The second boundary equation (3.62) is derived straightforwardly from a relation 
analogous to (3.52 ) using explicit r epres entation of T and V. Writing the scaling function 
Aiiyjj = A4(^5i; ^Ej) we expand (3.61) to 0(7i~ 2 ), in particular keeping the order 1/n 
coming from the anti-commutator of (3.50). Omitting straightforward details, we obtain 
a differential equation which fully determines A4(£i , £2), for £] < £ 2 

(d 2 + 7r^)y\4(fi,f 2 ) = ^ (P + (1 — 2£i)<9i + 2(1—^ 2 )^ 2 ) A4(^i, ^ 2 ), (3.63) 
M(0,£ 2 ) = A*(£ 2 ), <9^(0,£ 2 ) = 0. (3.64) 

This system is solved with an ansatz A4 (£i,£ 2 ) = C(£l,£ 2 ) + (^i)TW(^2) = 

cos(7t£i) cos(7r£ 2 )+C(£i,£ 2 )+n _1 cos(7r£!).M(£ 2 )+n _1 .M(£i) cos(7t£ 2 ) + (D(n~ 2 ) , resulting 
in an inhomogeneous system for the connected correlator C(£i, £2), with (3 cancelling out, 

(3.65) 

(3.66) 


71’ 


(d l +7T )C(£ 1,6) = —(£2 - l)cos(7r£ 1 )sin(7r£ 2 ), 


n 


C(0,£ 2 )=0, diC(0, £ 2 ) = 0, 
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.&) 


Figure 5. Scaled connected 2-point spin-spin correlation (3.67) function n x C(£i, £ 2 ) 
in NESS of XXX chain (A = 1). 


with a solution^ for & < £ 2 : C(fi,f 2 ) = “1^6(1 ~6) sin(7rfi) sin(7r£ 2 )- For & > 6, the 
solution is obtained from the symmetry C(£i,£ 2 ) — C(^ 2 ,^i), or generally (see Fig. [5]) 

7r 2 

F(£i,6) = «2^ mi n(6,6)(l - max(f 1 ,f 2 ))sin(7rfi)sin(7rf 2 ). (3.67) 

Note a qualitative resemblance to a 2-point function in classical SSEP (see e.g. Ref. [52]). 
apart from a trigonometric factor sin(-7r£i) sin(7r£)2) which seems to be of genuinely 
quantum nature. Our result establishes anti-correlation C < 0 between arbitrary pair 
of spins and the hydrodynamic scaling C oc 1 jn. Using this strategy one could derive 
further all the higher fc-point transverse spin correlation functions. 


3-4- Isotropic (XXX) chain - SU (2) — twisted boundary driving 


ffere we briefly discuss some details of explicit computation of a more general 
nonequilibrium partition function 

= ((0l| ® OT)T"(|tfR> ® |^)), (3.68) 


generalizing the expression (3.1), and observables in the XXX case with twisted 
boundary driving as described in subsect. 2.6.2 and treated originally in Ref. [62j. As 
using an arbitrary pair of twists requires full knowledge of the transfer operator T 
beyond the diagonal subspace /C we shall immediately utilize the rotational invariance, 
and choose a coordinate system in which the source axis is un-twisted while the sink 
axis is tilted in the x — z plane, specifically 0l = 4>r = $l = 0, = 0 6 [0,7r), where 6 

is the angle between the source/sink measurement axes. In such a case one can again 
work with the diagonal auxiliary subspace /C and the projected transfer operator T, 
rewriting Z n as (using asymptotic scaling of local transfer matrix elements (3.46)) 

2k 


Z n~J2 


k =0 


e 

tan - 
2 


T n \k) 


* Note a typo in the expression for C(£i,£ 2 ) in Ref- [70] . 
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iq(n !) 2 Zs ex P ( nf(k/n ) + 2/clog |tan0/2|) 


k =o 


— Fi(n\) 2i n / dt (dft/dt) exp {n(g t + 2^ t log |tan6*/2|)). 


(3.69) 


Note that k = 0 (y = 1) as required by solvability of twisted-driving boundary conditions 
(see subsect. 2.6.2[ ). This integral can be explicitly evaluated asymptotically (n —> oo) 
by means of the saddle point method, namely expanding around the extremum of the 
exponential at t — log tan( 0 / 2 ), yielding 


~ Fi(n !) 2 Sm ^ y/ 7 rn(l + ( 7 r - 9) cot 6) 

U 


2 n 


Z n -1 (tt-0 ) 2 


7T-y. \TI-0f ' v - A-* ■ ( 3 - 70 ) 

This reproduces the leading 1/e-order result of Ref. [03] and generalises it to any fixed 
value of e > 0 in the large n asymptotics. Note that Z n {6) is not continuous at 6 = 0, 


where one should take instead the expression (3.41), as the limits n —» oo and 6 —>■ 0 do 
not commute. 

Note that for the computation of Z n (9 ) described above we could still allow for some 
left-right driving asymmetry, and hence non-vanishing value of spectral parameter A, as 
long as x — 1- However, bellow we report, following Ref. [63], a simple calculation of 
vectorial spin-currents and spin-densities which is based on simple closed-form algebraic 
identities among T, V Q which are only possible for fully symmetric driving, i.e. A = 0 
and (3.49), so we assume this to be the case for the rest of this discussion. For 


SU (2)— symmetric XXX model, one can write the local conservation law for the full 
spin density vector a x = (ex*, crj, a 7 x ) and spin-current vector j x ,x+\ satisfying 


^ ®x [-H, a x ] j x — 1,2 


j x ,x+ 1 , where ] x>x+l := 2 a x x a x+1 , (3.71) 


and j x , x +i = \jxx +1 i s current discussed earlier. The expectation for the current 


components J a = (j xx+ f), which due to continuity equation (3.71) has to be site- 
independent, can be expressed in terms of the commutators 

r= ZlF) Z. WO, 01^, WIT”- 1 life,® (3.72) 

n /3,7S{x,y,z} 

where e Q/ 3 7 is Levi-Civita symbol. Facilitating algebraic identities, which hold for A = 0, 

[V x ,V y ] = 2i(t — s)T = 2iT(t — s), (3.73) 

[V y ,¥ z ] = (s + — s” + t + — t _ )T = T(s + — + t + — t~), (3.74) 

[V z , Y*] = (s+ + s" - t+ - t")T = T(s + + s~ - t+ -t~), (3.75) 


and elementary properties of the coherent states (2.103) we find the scalings of the 
in-plane current components 


r = 4i (t 


Z 9 9 

J x = 4i(g — t) n 1 tan - = — J z tan-, (3.76) 

Z n 2 2' 


n 2 

where ( = 4e/(e 2 + 4Zr). The transverse component behaves drastically differently 
though, it asymptotically scales in an ‘Ohmic’ fashion n -1 
2(t r-0) 


J y ~ 


n 


(3.77) 
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and behaves discontinuously at 9 = 0 where it vanishes (J y |gi =0 = 0). 

Vectorial spin-density profiles M x = (<j x ) are found in full analogy to computation 
described in subsect. 3.3, by extending the algebraic identity (3.50) to arbitrary 


components [T, [¥,¥“]] + 2{T, ¥“} = 4(s(s + 1) + t(t + 1))¥" due to SU( 2) invariance. 
Taking the continuum limit one again arrives to harmonic differential equation 
d 2 
,d £ 2 

for all three components of the continuous spin-density M a ((x — 1 )/{n — 1)) = M 
where appropriate boundary conditions follow from explicit representation of T, V" and 


—+ (vr-0) 2 AT(O = 0 


(3.78) 


a 
x ? 


the properties of coherent states (2.103), resulting in asymptotic harmonic profiles 


M z x — cos ( (tt — 6 ) 


x — 1 

I-7 

n — 1 


M x ~ sin ( ( 7 r — 9) 


x — 1 

)- 7 

n — 1 


Ml. ~ 0. 


(3.79) 


4. Hubbard chain 


Here we turn to a different, two-species quantum model, the fermionic Hubbard chain 
[23]. The Hubbard model is the fundamental model of strongly correlated electrons 
on regular lattices. Even though the model on a ID chain has been solved by the 
coordinate Bethe ansatz a while ago ra. it still poses many deep fundamental questions, 
in particular regarding its dynamical and nonequilibrium properties. Here we describe an 
explicit MPA solution of the corresponding nonequilibrium steady state of the Hubbard 
chain for diagonal (untwisted) boundary driving. We shall discuss graph theoretic 
interpretation of the solution and identify key elements of both approaches: IDO method 
(following Ref. [73|) and local operator divergence (or Lax operator) method (following 
Ref. [(75 ). 

The Hubbard Hamiltonian for an open chain of n sites, with canonical fcrmi 
operators c StX , x e (1... n}, s e (7, i}, reads 


H — 2 ^^( c l, x c s,x+i + c l, x +i c s,x) + u ^^(2n^ iX l)(2n± jX 1) 

S,X X 

+ AhXrat.i + n Li - !) + Ahtfat+ n i,n ~ 1), (4-1) 

where n S)X = c\ x c s x . The nondimensional interaction parameter u = U/(2t\f) contains 
standard Hubbard interaction U and hopping amplitude t^, while /xl/r are non- 
dimensional chemical potentials at the boundary sites which shall produce the coherent 
part of the boundary driving. The incoherent boundary driving is provided by four 
Lindblad channels which manifest a pure source/sink for electrons at rates Tl/r 


L\ — y/TYcji, Lo — V / Tl c {, 1 , -^3 — \ZTrC| )TI , L 4 — n/Trc^^. (4.2) 

We shall again be interested in the den sity operator p 0a of NESS defined as the solution 
of the stationary Lindblad equation (1.2), Cp^ = 0. Before proceeding, we shall 


reformulate the problem in terms of a spin-1/2 ladder, namely implementing the Wigner- 
Jordan transformation which expresses the anticommuting fermi variables 

r . — pD) fj- c , — p( CT ) p( r ) T ~ 

Wz — r x-l u x 1 — r n r x-l'x 


(4.3) 
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where Pj' 7 ' 1 := <J\<J 2 ■ ■ ■ cD, Pj r> := rfr| • • • r|, in terms of two sets of independent spins- 
1/2, a x , r x x E {1,..., n}, s,t E J = {+, —, 0, z}, a x = r° = 1, which can be considered 
as operators over / H® n . The local physical space is now four dimensional 7i v = C 2 ® C 2 
so that <j s T i span||]the complete basis of End('Hp). The Hubbard Hamiltonian (4.1) then 
maps to 

n —1 

H = h XtX+ 1 + /zl + hji, (4-4) 






hi,2 ■= K,2 + ^1,2 + g ( a l T l + a 2 T 2) , 

K ,2 '■= 2(7 1 °2 + 2<J 1 &2 > h l ,2 ■= 2r l T 2 + 2t 1 T 2 


U 


^L/R —O 


Z I ML/R / z | z \ 

1/nTl/n ' r\ y^l/n *" ^1/n) "> 


while the Lindblad jump operators map to 


n = Uny. u = Uvy'y, £ 3 = £4 = -v/Jv^pp-r. 


T )n-~ 


(4.5) 

(4.6) 

(4.7) 

(4.8) 


However, since the Hamiltonian and the dissipator T> = ^ x conserve the numbers 
of spin-up and spin-down electrons, N a = \ (o 7 x + 1), N T = Y^=i\i. T l + 4); 

[H, N a / T \ = 0, [N <t / t ,T>(p)\ = V([N a / T , p]), the unique steady state p^, should [TS] also 
conserve N c / T and their parities Pn^ T \ i.e., [p^, Pn^ T) ] = 0. Therefore, should also 
be a fixed point of £ = —iad H + Y^ a =i ^ L a , £poo = 0, where are replaced by 

L\ = x/Tlct 1 ”, £2 = \/rL r i + ) 4-3 = x/Tr^w , £4 = x/TRT n . (4.9) 

i.e., with all unitary conserved (non-local) operators removed (noting that (pO/O) 2 = 
1). This remarkable fact teaches us that non-locality of Wigner-Jordan transformation 
has no effect in the (nonequilibrium) steady state but potentially affects the nature of 
relaxation. Uniqueness of NESS can be again proved by straightforward application of 
the Evans-Frigeiro theorem [24, [28] trivially extending the argument of subsect, 
two species of spins. 


2.1 


to 


An important Z— symmetry of the Hubbard model, analogous to (2.19) for the 


Heisenberg model, is generated by the spin-flip operator G, i.e. permutation operator 
between cr and r spins (or fermion species), defined as Ga s G = r s , G 2 = 1. Clearly, 

Ghl 2 G — h\ 2 , Gh h2 G = hi, 2 , GHG = H, GV(p)G = V{GpG). (4.10) 


4-1. Walking graph state representation of NESS 

In the absence of previously known non-Hermitian Lax operators with enough free 
complex parameters for the Hubbard model (note that the Hermitian Shastry’s Lax 
matrix [83] would not work as it lacks a free representation parameter), we shall 
again start with a constructive approach of IDO method [73], while impatient reader is 


welcome to jump right away to a more elegant formulation of subsect. 4.2 


jj Note that here we use letters s, t to name indices denoting physical space components in 
contradistinction to previous sections where they denoted complex spin parameters. 
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A useful technical result which can be implemented to establish NESS is a trivial 
extension of Lemma 1 to a symmetrically boundary driven Hubbard ladder: 


Lemma 2. Let £ End(%p n ) satisfy the following conditions (defining relations): 
(i) a recursion identity for the bulk, setting Hi := I 4 , 


[H,Q n \ = 


-is 


(a z r s 0 I 5 , + u 0 P s n L u ! - 


s,0 


yO,—S . 


se{o,+} 


)a r 


Qn S - 1 °' 


)a~ s r z 


(4.11) 


introducing the operators Pf'lj , Q'lfii 6 End('Hp (,t 


p s ,t _ 

r n -1 — 


tr 1 {(c r i'7'i) t H ra } 


Qn-1 = 


tr n{(<#7-*) t n n } 


tr({<r s r t )i(T s r*} ’ 1 tr{(er s T t )'lcr s T t } 

and (ii) the boundary conditions (rendering fl n upper-triangular with unit-diagonal) 


(4.12) 


Pffl 1 = 0 for s £ {-, z} or t £ {-, z}, 

Qn-i = 0 for s e {+, z} or t £ {+, z}. (4.13) 


Then, the density operator 

Poo = Roc = , (4.14) 

satisfies the fixed point (NESS) condition 

4 

i[#,Poc] = ^PlAPoo) (4.15) 

a=l 

for symmetric, totally incoherent driving £ = T L = r R , Pl = Pr = 0. 


A straightforward proof along extension of Eqs. (2.12 2.17) is left to the reader. 

The following constructive strategy for obtaining exact NESS solution has been 
devised [73] which is based purely on empirical data about the model. One starts by 
computing numerical NESS density operators for small systems, feasible for n < 6, 
and determine their Cholesky factors Q n . Operators Q n posses U( 1) x 77(1) symmetry, 
namely they commute separately with the species number operators [Vt n ,N a / T ] = 0 
hence all non-vanishing terms of a general operator expansion Q n = t c s,t &"=! cr Sx T tx 
should satisfy Y^=i d(s x ) = 0 and d{t x ) = 0 where the shift-function d : J —» Z 

is defined as d(±) = ±l,d(0) = d( z) = 0. Thus each sequence (si,H,... ,s n ,t n ) e J 2n 
with c s t 0 can be considered as an n —step recurrent walk on a 2-dimensional cartesian 


grid Z x Z originating from site (0,0), visiting a point YT y =i(d{sy), d(t y )) after step x. 
However, empirical evidence suggests that the set of non-vanishing terms is much more 
restricted and can be compactly encoded by a directed graph (V, £) depicted in Fig. [ 6 j 
The set of all visited nodes (or vertices) V C Z x Z is composed of: the origin 
0 = ( 0 , 0 ), the diagonal nodes k = (k,k), and upper-, and lower-diagonal nodes, 
(k — 7)+ = (h — 1 , k), and (k — 7)~ = (p, k — 1), for k G {1, 2,...}. Cartesian coordinates 
of a node v = (rd,u 2 ) will be written as v u , v £ {1,2}, in general. The set of directed 
edges S(G) contains: vertical, horizontal, diagonal, skew-diagonal, and self-connections, 
as indicated in Fig. [ 6 j where only self-connections of diagonal nodes are degenerate 








Matrix product solutions of boundary driven quantum chains 


39 



Figure 6. (from 73]) Diagram of a semi-infinite graph G (structure repeating 
periodically beyond the upper-right corner) showing the allowed transitions for building 
up the MPA form of NESS for the Hubbard chain. Nodes in black, edges with 
multiplicity 1 in red, and edges with multiplicity 2 in blue. Each edge e is associated 
with a physical product-operator w(e) = a b r b ~ where b 1 ' = 0 ( b v = z) for edges 
connecting white (black) nodes, where v is the Cartesian component which does not 
change along such e in the diagram. Degenerate edges correspond to operators a °r° 
(y = +1) and ct z t z (/z = —1). Insets indicate all possible terms (two in each, orange 
and brown) for two examples of [h,w(e) <B> <*;(/)], specifically [h,a + T + <g> cr + r°] (a), 
and [h,a°r~ <B> ct 0 t 0 ] (b). Full arrows denote valid edge factors, while dashed arrows 
correspond to defect operators. 


with multiplicity two. Edges may also be identified with triples e = (p(e), q(e); //(e)), 
pointing from node p(e) to node q(e) and having degeneracy label /z(e), where pt = 1 for 
all edges except diagonal self-connections (. k,k;n ) where pt e {±1}. 

There are two crucial concepts of the IDO method generalising the concept of 
defect operator a z in the XXZ model. The first is the index function cu : £ -A End('Hp) 
associating a local physical operator oj{e) to each edge e of the graph, which can be 
fully determined by a careful inspection of empirical data, i.e., it should match a Sx r tx 
for the edge corresponding to x —th step of all walks generated by nonvanishing terms 
of Q, n . Painting the nodes of the graph as black, white and black&white (see Fig. [6j) we 
encode the empirical data suggesting the index function 

w(e) = c/ (e V 62(e) , where b u :S^J, (4.16) 

as follows: b u (e) — ± if q u (e) — p v (e) = ±1, while for q u (e) = p u (e), Ifle) = 0, 
if e connects white nodes, and b u (e) = z, if e connects black nodes. For diagonal 
self-connections (on black&white nodes), the index functions are determined by the 
degeneracy index, b u (k : k; 1) = 0, b u (k , k\ —1) = z. 
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The second key concept is the defect edge of the graph. Let us consider an arbitrary 
walk of length 2, i.e., a pair of subsequent edges e, / G £, with q{e) = p(f)- Writing a 
Hamiltonian density on a pair of sites as hi j2 = hf 2 + ^12 + u i a i r i + U 2 ^ 2 t 2 , which can 
represent either the bulk or boundary part of H (4.4), one finds the following general 
form of the local commutator of h with a tensor product of two valid edge factors for a 
pair of consecutive edges (2—walks) e, / G S, q(e) = p(f) 


[h,uj(e)®oj(f)} = 


p(e')=p{e),q{f)-q(e')=d(s,t) q(f')=q(f),p( e )-P(f)= d ( s ’ t ) 

s,tej,e'e£ s,teJJ'e£ 

r S,t 1 


(4.17) 


with suitable structure constants X^j(ui,u 2 ),Yfj(ui,u 2 ). We define a displacement 
vector associated with a pair of Pauli indices as d(s,t ) = ( d(s),d(t )). Eq. (4.17), in 
analogy to identity (2.23) for XXZ model, has the following crucial property: Any 
tensor factor cHr* in the first (or second) sum on RHS of (4.17) is (i) neither of the 
form u ;(/') (or uj(e')), for any edge f (or e 7 ) which would complete the 2-walk (e 7 , /') 
to connect the same nodes as (e, /), (ii) nor is the missing link d(s,t) between q(e') and 
q(f) (or p(e) and p(f')) provided by any edge of the graph. We shall call such a factor 
a defect operator, or defect edge if referring to the graph. See insets of Fig. [6] for two 
characteristic examples. 

To each vertex v G V we associate a vector space H v , such that the entire auxiliary 
space is a direct sum "H a = ©,. eV l~i v , and associate a transition amplitude to each edge 
e G £ as a linear operator a e G Lin ('H p (e),'H q (e))- Writing a WGS ansatz (2.37) for the 
amplitude operator of NESS 


- - n 


^ a ei a £2 ■ ■ ■ a £n uj(ei) ® u(e 2 ) 0 ■ ■ ■ u{e n ). 

(ei,..,e„)6W„(0,0) 


(4.18) 


one notes that, since the Hamiltonian (4.4) is a sum of local terms, the entire commutator 
[H, Q n ] written in the tensor product expansion (like (4.18)) is composed of terms which 
correspond to n-walks over a defective graph with exactly one defect operator. As 
the RHS of (4.11) has only boundary defects, in the first or last factor, all the terms 
with defects in the bulk should therefore identically vanish. Picking any pair of nodes, 
v, r G V, which can be connected with at least one 3—walk, it is then sufficient that the 
following local conditions are satisfied 

[H 3 ,u(e) ® u(f) 


cole 


a t 


i u 


= 0, (4.19) 


( e,f,g)eW 3 (v,r) 

for any pair of edges e',g' G S(G) for which p(e') = v,q(g') = r, and any defect 
component s,t G J. Here and below G End(7^p ) denotes the Hamiltonian (4.4) 
for a small cluster of n = k sites with open boundaries Pl/r = 0. Of course, for 
many combinations (v,r,e ', g ', s,t) the above equation is trivial, i.e. always satisfied, 
e.g., when cr s r 7 = tv(f') for some valid edge /' between q{e') and p(g'). The remaining 
local equations which need to be satisfied to yield ( 4.11[ ) are those for which the defect 
operator sits at the first x = 1 or the last x = n tensor factor. Again, one can factor 
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out sufficient local conditions, which can now be formulated on two sites, in terms of 
2—walks, namely 

y a e a/tr| (cr s k 0 a;(/ / )) t (\H 2 ,u(e) ®w(/)] +i eP{u(e)) ®<u(/))} = 0, (4.20) 
(e,/)ew 2 ( 0,«) 

a e a/tr| (co(e') 

(e,f)£W 2 (v,0) 


—S—t\ t 

a r 


) f ([H 2 ,u(e) ® u(f)] -iau(e) <B> P(a>(/))) } = 0, (4.21) 


for all e', f G £, with q(f) = v, and p(e') = v. V is a, map over End('Hp) defined as 
V(p) := \cr z <8)tr CT (p) + |tr T (p) where tr CT (or tr T ) denotes the partial trace over cr (or 
r) qubit of LL P . Here, the set of possible defect operators is quite limited, specifically, 
to (s,t) € {(0, z), (z, 0), (+, z), (z,+)} for the left boundary conditions (4.20), or to 
(s,t) e {(0, z), (z, 0), (—, z), (z, —)} for the right boundary condition (4.21). 


Summarizing, finding a e obeying the three-point recurrences in the bulk (4.19) 
with the two-point boundary conditions ( 4.20[4.21 ) is sufficient for establishing validity 
of Eq. (4.11) with ansatz (4.18) together with the conjectured structure of the graph 
(V, S) and its index function cu and hence, according to Lemma 2, exactly solving NESS 
for symmetric incoherent driving Tl = Tr = e, /zl/r = 0, for anyn. The solution, unique 
up to gauge transformations, has been found [73] by means of a computer program in 
Mathematica, requiring the local auxiliary spaces of uni-color nodes 0, {k — I) 1 * 1 to be 
one-dimensional and those of black&white nodes k, k > 1, to be two-dimensional: 

O(0,0;+l) = 1) O(0,0;-l) = 0, O(0,l/2±) = — G ®(l/2±,0) = — h 

\k i 


a ((fc-l/2)±,(fc~l/2)±) — —(-1) 

O(o,i) — ( —2ie 0), a(i,o) - 


a ((fc-l/2)±,(fc-l/2)T) — 1£, 

1 


2 e u 


a (k,(k+l/2)±) — 


0 


a (k,(k- 1/2)±) — 


-(-l) k (ku + k)f 


<*((*—i/ 2 )±,k) = ( ~£ 0 ), a^ k+1/2 )±,k) = (—(—1) (i + fe«) (-l) L2j f), 


-1)L^ J | 

k <■- i ^ /'_iUVds 


a (fc,fc+i) — 
a (fe+i,fe) = 


-l) fc 2ie 


1 0 
0 0 


a 




l (fc,fc;-(-l) fc ) 


-l) fc (i + %eu)(i(k+l)u 
(-l)^j (±keu - 1) 
(-l) fc (l - \keu) 

0 


-1) [ 


^ 1J (| - i(k+l)u% 


-l)L^Jk 


¥ 


(-!)' (| -iM| 


- 1 ) 


kJi £ 


(4.22) 


4-2. Lax representation of NESS 

Having an exact NESS solution for the Hubbard model at hand, one can now 
explore its mathematical properties more deeply. A strong motivation for that 
comes from observation that MPA formulation of the WGS ansatz ( 4.18[4.22 ) Q n = 
t ■ ■ ■ A s n ,t n cr rk 1 • • • a n T n n with A s,t = ® ee£: ^ s.b'(eft, b ^ e) a e allows an explicit 
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factorization A Sjt = S s T*X, with S s .T,,X G End(H a ), [S s ,T f ] = 0, which is, 
in spirit, very close to Shastry’s form of the Lax operator [83]. We show in this 
section how a general Lax form of the NESS amplitude operator Q n can be derived 
satisfying a generalised Sutherland relation (closely following Ref. [65] ) which reproduces 
and generalizes the result of the previous section, namely it solves the boundary 
driven ffubbard chain for arbitrary rates T L / R and chemical potentials p l/r. IDO 
technique could then be re-interpreted merely as a graph theoretical representation 
of the Sutherland condition formulated locally between adjacent vertices of the graph. 

It turns advantageous here to choose a particular basis of auxiliary sub-spaces for 
diagonal (black&white) vertices TLk>\ = lsp{|fc _ ), |fc + )}, 'Ho = lsp{|0 + )}, and hence 
to identify the nodes of the graph with unique labels of individual auxiliary basis 
states V = {0 + , ^ + , | , 1 _ , 1 + , | + , § , 2~, 2 + . . .}, so that the entire infinite-dimensional 
auxiliary space is a simple linear span H a = lsp{|r>) ]v G V}. We extend the definition 
of the spin-flip G over H a as a diagonal reflection of the graph, G | A: ± ) = |/c ± ) , 
G \k+^ ± ) = |&; + | T ), k G Z + . We begin our analysis with a simple observation: 

Lemma 3. TUBf Assume there exist operators S,S,S,T,T,T G End(H a <S> H p ), and 
X, Y G End(H a ) (acting as scalars over Up), satisfying 


[hi 2 , SiXSa] = S!XS 2 - SiXSa, (4.23) 

[hf 2 . TiXT 2 ] = T x XT 2 - ThXT,, (4.24) 

ST + TS — ST — TS = [Y — ua z r z , ST], (4.25) 

[S, T] = 0, (4.26) 

[X, Y] = 0. (4.27) 


Subscripts, like in S x , indicate independent local physical spaces pertaining to sites x in 
the embeded representation End(H a <g> H p n ). Then, one can define a Lax operator and 
its ‘derivative’ L.L G End(H a ® H p ) as 


L = STX, (4.28) 

L = 1(ST + TS + ST + TS — {Y, ST})X, (4.29) 

such that the following Sutherland-Shastry relation (or generalized local operator 
divergence condition) holds 

[hp 2 , LjLa] = (Li + YL0L 2 - LRL, + L 2 Y). (4.30) 


The proof is a straightforward insertion of ( 4.28[4.29[ ) into Eq. (4.30) followed by 
subsequent application of identities (4.23 4.27) observing the definition (4.5). 

We continue by deriving an explicit closed form representation of algebraic identities 
(4.23 4.27). Assuming the spin-flip symmetry 

GSG = T, GSG = T, GSG = T, [G, X] = [G, Y] = 0, (4.31) 


and writing out the components S = 
S,S.T.T, we find that Eqs. 


s£j_ 


S s a s , 


T = Ylte 7 T f H, and similarly for 


(4.23) and (4.24) are equivalent, Eq. (4.25) is invariant 
under G, while Eq. (4.26) implies [S s . T 4 ] = 0. Eqs. ( 4.23[4.24 ) are in fact just 
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a particularly ‘decorated’ 6-vertex Yang-Baxter equations for free fermion (or XX) 
chains. We shall thus make an ansatz for Sb T f in which each square plaquette 
{k + , k +| + , k+\ , k+l~} of the graph spans a pair of representations of a free fermion 
algebra (see Fig. [7]), namely requiring that {S + ,S“} (and similarly for {T + ,T - } via 
(4.31)) is in the center of the algebra generated by S s , T f 

[{S + , S~}, S s ] = [{S + , S-}, T f ] = 0, s,tej. (4.32) 

One finds that these conditions are fulfilled by an ansatz 


S + = \/2 (j^ + ) (k + l ' I + 1^+1 ) (k + 1 ), 

k =o 

oo 

S = \/2^^(—l) fc [\k-\-\ ) (k + \ + \k+l ) (k+\ , 

k =o 

OO 

S° = ]T(|2fc+) (2k + \ + |2 k+\ + ) (: 2k+\ + \ 

k =0 

+ \2k+l ) (2/c+l | + |2/c+^ ) (2fc+^ |) 

oo 

+ (|2fc-i + ) <2fc-fl + |2 k-) (2k-\) , 

k =1 
oo 

S z = ^(|2£;-1 + ) (2k—l + \ + |2 k-\ + ) (2k-\ + \ 


(4.33) 


k =1 


+ \2k~) (2k-\ + \2k + l )(2 k+\ |) 


oo 

+ A^ (\2k + \ + ) (:2 k + \ + \ + |2fc+l") (2A:+1“|) , 


k =0 


where A G C is a free parameter. Eqs. (4.31) imply definition of another set of auxilliary 
fermi operators = GS 4 G, such that Eq. (4.26) is satisfied. 

Furthermore, one can write a consistent ansatz for the ‘interaction’ operator X 
coupling the neighbouring plaquettes: 

OO 

x = |o+) (o+i+£(-i)* 5^ imr'Yi 

vye{-,+} 


k =1 


oo 

+ w ^(—1) A (\k+\ ) {k + l l + l k + \ ) (k + l I j , 


(4.34) 


k =0 


where X k = {Xf’ u } I/ y e f >+ \ are still unknown 2x2 matrices and w E C is another 
free parameter. Namely, Eq. (4.23) yields a system of linear equations for auxiliary 

(4.35) 


operators S S X, XS S , with a unique solution parametrised by X k , w, A: 

OO 

S+X = -2v / 2^(-l)W+- I*-) (fc+pl, 


k =1 
oo 


S-X=-2v / 2j].Y t - + |fc + ){fc-y|, 


k =1 
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XS+ = 2v / 2]T(-1 t x t \k-D (4+1 

k =1 
oo 

XS- = -2V2j2 X k + \k+\ + ) (k~ I, 

k =1 
oo 

S°X = XS° = 2^(w|2A;-l + ) (2fc-l + | - w\2k~) (2 k~ 


k= 1 
l + \ 


- |2fe-D (24-§+| - |24-| ) (24—| |) 

OO 

2A^(-«> 24+) (24+| + XJ,;, |2*+1“> <2fc+i'|). 


k=0 


S Z X = XS Z = 2 |2fc+l“) (2fc+l‘| - w \2k+) (2 


k =0 

l+\ 


+ +^2fc + l 2 ^+| )( 2 ^+^ l+^2fc+ll 2 ^+i )( 2 ^+| I) 

OO 

+ 2A^(«,|24-l+> (24-1+1 - X 2 y |24-n <24-p|). 


k =1 


Assuming X to be invertible (i.e., w^O, clet X k ^ 0) and plugging expressions (4.35) to 
the remaining identity (4.25) result in (i) a unique consistent expression for the ‘spectral’ 
operator Y 


Y = -2\u^y\k + ) (k + 1 + \k + l~) {k+l~\), 


(4.36) 


fc =0 


which clearly commutes with X, as required by (4.27), and (ii) recurrence relations 
for the matrix elements of X k \ Xf^ = Xft~ — uw, Xf+ x = Xf + — uw( 1 — A 2 ), and 
det X k = —w 2 , while also fixing the initial condition X^ + = 1, X Q = —w 2 , yielding 


X k (X,w) = 


(w + ku)w 
—kuw 


(4.37) 


1 — (w + ku)w( 1 — A 2 ) 

1 — kuw (l — A 2 ) 

Note that Xf + /Xf~ can be chosen freely exploring a gauge freedom jAA) —* ^ ±1 | ), 

k = 1,2... We have thus constructed two-parameter representation of the Lax matrix 
L(A,w) = S(A)T(A)X(A, w) satisfying Sutherland-Shastry relation (4.30). We propose 
to call A a spectral parameter and w a representation parameter. Remarkably, our 
representation is generically of infinite-dimension, for any nonzero u, and seems to be 
essentially different from Shastry’s [83] which, including appropriate auxiliary operators 
S, S,..., forms a 4-dimensional representation of the algebra (4.23 4.27). 

As an application of such novel Lax operator we consider a markovian master 
equation and demonstrate how to obtain NESS fixed point (4.15) for general 
Hamiltonian (4.4) and boundary dissipation (4.9): 

Theorem 1. 165J Unique fixed point Cp^ = 0 of boundary driven Hubbard chain reads 
Poo = (trR 00 )~ 1 R oc , Roc = tt&K, (4.38) 
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Figure 7. Diagrammatic representation of factors of the Lax operator where 
auxiliary states are labelled by vertices V. Diagrams for T s are obtained by 
reflection of those of S s across the diagonal. Red/blue arrows indicate offdiagonal 
transitions with amplitude ±\/2. Red, blue, green, black, open points represent 
diagonal multiplications by w, —w, oc A, 1, 0, respectively, and brown circles represent 
multiplications by 2 x 2 matrices X^. 


where = Q n (A,w) is a highest-weight transfer matrix 
n = (0+1 L X (A, iu)L 2 (A, w) ■ ■ -L n (A,w) |0+) 
and Ii is a diagonal operator 

K = I\] K, ■ ■ • K n , K n = exp + r|)) 


(4.39) 


(4.40) 


with k = 4 logr L /r R and parameters A ,w are related to coherent and incoherent biases 

^ - r - r —-c, w — - (^l - +1 (r L + r R )). (4.4i) 

r L + r R -i(/t L -hr) 4 

Proof. Let us now invoke two copies of the auxiliary space and define operators 
S, T. S'. T' e End('H a <g> ® H p ) as 

S = S s ® l a <g) cr s , T = T' <8) l a <g) t 4 , and 


s'=E^ 




0 T , T' = ^l a 


(r 4 ) 


t\T 


() T denotes the matrix transposition and S the complex conjugation, i.e. replacement 
A, ta, —>■ A,w, and similarly for S, S, S', S', T. T'. T', and X, X',Y, Y' 6 End("H a <g> 
"H a ). In fact, the primed operators S', S', S', T'. T', T', X', Y' generate a conjugate 
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representation of the algebra (4.23 4.27). Noting [hi i2 , K\K-f\ = 0 and the Jacobi identity 
one finds that the following double auxiliary operators 

L j = LjL'jKj, L j = (LjL'j - I. L')/\,. Y = Y - Y', (4.42) 


also respect Sutherland-Shastri relation (4.30), resulting in the telescoping series 

n— 1 

[hjj+i, LiL 2 • • • L n ] = (Li + {Y, Li})L 2 • • • L n 

3 = 1 

— Li • ■ •L n _i(L, t + {Y, L n }). 

Double Lax operator expresses NESS in a compact form 

Roo = (0 + , O+IL.L, • • • L,„ |0+, 0+), 


(4.43) 


(4.44) 


hence the fixed point condition CR = 0 becomes, after applying ( 4.43[ ) to [H,R], 
equivalent to a pair of equations for ultralocal operators at the boundary physical sites 


(4.45) 


(0 + , 0 + 1 (ir L (V a+ + 2} r +)IL + L + LY + [hi,, 1L] J — 0, 
ir R (P CT - + V T - )L - L - YL + [h R , L]) |0 + , 0 + ) = 0, 


where boundary interactions with fields, /il/r, are defined in (4.7). Using explicit forms 


(4.334.37) and in particular X|0 + ,0 + ) = X'|0 + ,0 + ) = |0 + ,0 + ), each of Eqs. (4.45) 
results in dim R p x dim R p = 16 equations for (bra/ket) vectors from R a ® 77 a , most 
of them trivially satisfied, whereas the non-trivial ones are equivalent to conditions 

dOlb. n 


We have presented infinite-dimensional irreducible representation of Lax operator 
and Sutherland-Shastry compatibility condition and shown how it can be employed to 
yield exact NESS of asymmetrically boundary driven Hubbard chain with arbitrary 


boundary chemical potentials. We shall discuss later in subsection |6.4| how this Lax 
operator can be used to define new conserved operators of the Hubbard model which 
break particle-hole or spin-reversal symmetry, in full analogy with the situation in the 
XXZ model. One could now also embark on computation of local observables in the 


Hubbard model. For example, linear dependence of the amplitudes (4.37) on auxiliary 


state k immediately yields, following the same asymptotic analysis as for computing the 


nonequilibrium partition function in the case of XXX chain presented in subsect. 3.2 
a universal scaling of the spin/charge currents J ~ n~ 2 and cosine-shaped spin/charge 
density profile, as observed in numerical simulations of Ref. 


5. Lai-Sutherland spin-1 chain 

As the third characteristic example we outline (following Ref. [36]) exact solution of a 
boundary driven spin-1 Lai-Sutherland chain which, due to the three-state nature of the 
model with only a pair of states being dissipatively transformed at the boundary, allows 
for macroscopically degenerate NESS. 
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Figure 8. Schematic illustration of the degenerate boundary driven 3-color quantum 
chain, or spin-1 Lai-Sutherland chain. The incoherent boundary jump processes 
transform just a pair of colour states (red and blue), while the green states remain 
invariant and hence the number of green particles is a constant of motion. 


Consider a finite chain of n sites with 3—state local physical space "H p = C 3 . Using 
the Weyl matrix basis {e u = \i)(j\ ; ,i,j = 1, 2, 3} of End ('H p ) = 0l 3 , we define a full set 
of local generators of the full matrix algebra £ = End('H n ), where = "H p n denotes 
3 n -dimensional Hilbert space of n-site chain, as 

(5.1) 


satisfying the Lie algebra relations 


= (S jh e“ - 6^)6,, 


i l 




(5.2) 


The spin-1 Lai-Sutherland model [50] for a chain of n sites is given by the Hamiltonian 

n— 1 

H ^ ^ h Xj x- 1-1) h X: x- (-1 $x ' d'/: - 1 T ($x ‘ Sa;+l) 1) ('-*•3) 

x=l 


where s x = 1 ^ ° 2 ° 3 


s x> s x> s xh with 


s 1 = —(e 12 +e 21 +e 23 +e 32 ) s 2 = — 

°x y/2' X ’ Cx ’ Cx ' Cx b ° x ^2 


(e 21 —ei, 2 +e 32 —e 23 ) 


si = el 1 - e; 


.33 


(5.4) 


form independent spin-1 variables (local s — 1 representations of SU 2 ) satisfying 


[ S a;) S x'] ~~ * '/ > y e ijk s x fix,x' ■ (5-5) 

k 

Straightforward inspection shows that the local Hamiltonian h X)X+ \ - the interaction - 
is in fact just the permutation operator between neighbouring sites 

hx,x+ 1 = 1 3 0r ” 1) ® I i,j) CmI ® lf (n_x_1) = e* J 4+i- ( 5 - 6 ) 

i,j =1 *,3=1 


The local Hilbert state basis is therefore given by a triple of states |1) = |t), |2) = 
|0), |3) = |J,), which can be interpreted as three different particle components (or 
colours); respectively, as spin-up (red) particles, spin zero or holes (green), and spin- 
down (blue) particles. 

Since Lai-Sutherland chain is a multi-colour quantum model one may associate 
with it a skew-symmetric tensor of particle currents, with two-site density 

J ij = i{e ij ®e ji -e ji ®e ij ), J l J = 1®°* _1) ® (5.7) 
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which, by construction, satisfies the following continuity equation 


d 

df 




ip, 4 


pJJ] — po 
'-'X J U X—1,X 


TO 

° x,x+l ' 


(5.8) 


can be considered as a partial current of the particles of colour i into particles of 
colour j. The total current of component (colour) i, 

3 

j = E ( 5 - 9 > 

3 = 1 

then fulfills the continuity equation 



P — V 

° X— 1,X ' J x,x+ll 


(5.10) 


where e“ can be considered as the operator of particle density of colour i. 

We shall now open the Lai-Sutherland chain and couple it to the environment via 
Markovian processes which act only on local quantum spin spaces at the boundary, via 
the following Lindblad jump operators 


L 1 = \leeY’ 



y(«J 2 , where s± := 4±i4- ( 5 - n ) 


Two dissipation channels, interpreted as the left and right magnetization bath, perform 
the incoherent processes |j") —* |J,) and |j,) —>■ |t), respectively, with the rates e. Both 
processes keep the hole state |0) unaffected. Since also the bulk dynamics generated by 
Cq conserves the number of particles of each colour, it follows that the whole Liouvillian 
dynamics (master equation) preserves the number of holes. More precisely, defining the 
hole-number operator N 0 e ^ as 


No |* 1 , * 2 , ■ ■ ■ , in) 



|* 1 , * 2 ) ■ ■ ■ ; in) i 


(5.12) 


we have that the set of all, Hamiltonian and jump operators, commute with No 


[H, N 0 ]— 0, [T 1)2 , N 0 ) — 0, (5.13) 

which implies that N 0 generates a strong [12] U( 1) symmetry of the Liouvillian flow 
(0. Nq foliates the physical Hilbert space into n + 1 orthogonal IVo-eigenspaces, 


Tin = 0" =o TLn \ NoHn 1 = oHlf 1 ■ The theorem A.l of Ref. [T2j then guarantees that 
the full Lindblad dynamics (1.2) is closed on — End ('Hn ^), and that 

a fixed point p ^ = lim^oo exp(t£^)pQ 1 '- 1 (NESS) exists for each symmetry subspace 
flow 


/H — 


H 


= -i + (N + T>l 2 )p 1 P = o. (5.14) 

The theorem by Evans and Frigeiro [241 [2B] can then again be used to show uniqueness of 
NESS p { £ for each fixed v. In the following we shall outline a simple algebraic procedure 
for actual explicit construction of density operators p&. 






49 


Matrix product solutions of boundary driven quantum chains 


5.1. Degenerate matrix product solution 


Let G End (J) be a projector to & v \ orthogonal with respect to Hilbert-Schmidt 
inner product {A\B) = tr A^B, with respect to which the Weyl basis e* ui 0 e l2T2 ■ ■ ■ is 
orthonormal. We define a grand density matrix of NESS as a direct sum of non-trivial 
solutions of (5. 14[) for all u, 


o = V o {u) 

poo / Poo 1 
v=0 


with pr = p <, 'Voo f o, 


(5.15) 


being solution of the fixed point equation (5.14) as well. The grand state p^ shall be 
sought for in terms of Cholesky factorization (in analogy to previous solutions of XXZ 
and Hubbard models) 

Poo(^) = ^n(e)^n(e), (5.16) 


where O n (e) G End (H n ) is some yet unknown operator which is represented by an 
upper triangular matrix in the computational basis \i \,..., i n ). Introducing an auxiliary 
Hilbert space - separable, but of infinite-dimensionality as will become clear later 
- we define the monodromy operator M(e) G End ("H„ 0 "H a ) as a spatially-ordered 
product of some local Lax operators L x (e) G End (fH n 0 "Ha), 


M(e) = L 1 (e)L 2 (e)---L n (e). (5.17) 

Index free Lax operator is defined as L(e) G End (Tip 0 B a ) so that one writes 
L x (e) = 1®0 L(e) 0 Furthermore, we define the components of Lax 

matrix L u "(e) G End('H a ), such that 

3 3 

L*(e) = Y 4 ® L*(e), Me) = Y e” 0 L*(e). (5.18) 

i,j =1 i,j =1 

We further assume existence of a special state |vac) G "H a , such that Cholesky factor 
writes as the auxiliary expectation value of monodromy operator, or equivalently, as 
MPA 


Q n = (vac| M |vac) = Y^ (vac| L* Ul • ■ • L*"- 7 " |vac) e* Ul 0 ■ ■ ■ 0 e*"- 7 ". (5.19) 

Fixing an arbitrary, fixed orthonormal basis {| i/>k)} of B a we define the conjugate 
Lax matrices L(e) by (^|L y (e)| ipi) := (ipk\ L lJ ’(e) \ifi). For notational convenience 
we denote the second copy of auxiliary space carrying conjugate representation of L 7 
as H a . One can then write MPA for NESS density operator directly, by introducing 
two-leg Lax matrices L u (e) G End (fH a 0 ?f a ), and L x (e) G End (fH n 0 B a 0 H a ) as 

V j (e) = Y M k {e) 0 L j \e), L x (e) = Y 4 ® L ij ’(e), (5-20) 

k i,j 


namely 


(5.21) 
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Note the transposition in the quantum space of the conjugated factor of (5.20). Here a 
two-leg monodromy operator 

M(e) = L 1 (e) • • • L„(e) e End (H n 0 0 TL a ), (5.22) 

and a product of a pair of vacua ((vac| = (vac| 0 (vac|, |vac)) = |vac) 0 |vac) have 
been introduced, so that (5.21) is merely a formal rewriting of (5.16). These definitions 
become particularly handy when we consider evaluation of expectation values of local 
observables with respect to NESS Poo{s). 

Let rj : = ie be a complex-rotated coupling parameter and let us (for convenience) 
relabel the quantum space matrix elements of the L-operator as 

l f t + v H 

1° u+ ] . (5.23) 



u 


u 1 

P 


Explicit MPA structure of the grand density operator of NESS is then established by 
the following result [36] : 

Theorem 2. 1361 (i) Suppose that 9 matrix elements {L' J } generate the Lie algebra g 
defined by commutation relations, 

[u + ,t ± ] = [u-,t ± ] = [u ± ,v ± ] = [t ± ,v ± ] = 0, 

r,u ± ] = [i^t ± ] = [iM^] = o, 

[1^] = ^, [P,^] = TVU k , 

[u =t ,v T ] = ±?/t T , [t =b ,v T ] = 
r, v ± ] = [P, v*] = =F?/v ± , [v+, V"] = 7?(P 

[t + ,t~] = [u+,u“] = 77I 0 , 

[1^,1°] = [ 0 =*=, 1°] = [v ± , 1°] = [t ± , 1°] = 0, 

with a representation over the Hilbert space H a satisfying the following conditions 
l 1 |vac) = 1° |vac) = l 1 ' |vac) = |vac), 


rj, 


(5.24) 


(vac 


P = (vac 1 1° = (vac 1 P = 


,vac 


t + | vac) = u + | vac) = v + |vac) = 0 , 

(vac1t = (vac| u~ = (vac| = 0. (5.25) 

Then, the grand state solution (5.15) to NESS fixed point condition (5.14) is given via 
Cholesky factorization (5.16) with explicit MPA (5.19) for f2 n (e) with rj — ie. 

(ii) A possible irreducible explicit representation of Lie algebra 53 (5.24) satisfying 


(5.25) is given as 

t + = b 




- .mUt 

= b j 


t = r/bl, 


4 -’ 


= r)b h 

v + = ?/(b t b| + s + ), v” = ?/(bjbj 

1 U = 9 (b{^b u + ^ - s z ^ , 1° = 1, (5.26) 

in terms of three auxiliary degrees of freedom with a three dimensional lattice 
{| j,k,l ), j,k,l G Z + } forming a basis ofTL a , namely, two bosonic modes b^ 

b{ |i, k, l ) = yjj + 1 |j + 1, k, l ) , b t |i, k,l) = y/j |i - 1, k, l) , 
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b l b\M) = Vk+ 1 \j,k + 1,1 ), \j, k,l) = Vk\j, k-l,l), (5.27) 

and a complex spin (Verma module of sl 2 ) 

s + | j, k,l) =1 1 j, k,l- 1), 
s“ | j, k, l) = (2 p - l) | j, k, l + 1), 

s z | j,k,l) =(p-l)\j,k,l). (5.28) 


with |vac) = |0, 0, 0) being the highest-weight-state. The complex spin parameter p should 
be linked to dissipation parameter via 


V = 


1 

2 


1 

V 



(5.29) 


Proof. The proof is based on verifying that the Lie algebra g, given by (5.24), can be 
equivalently defined by means of an appropriate Sutherland relation 


\hx,x + ll ^x^x+l ] = B x L, +1 - L X B X+1 , (5.30) 

with the-so-called boundary operator B x {e) e End {TL n <8> Ti & ) - operating non-trivially 
only in the local quantum space 


B x = q (< 


,33 


la - e 


11 


la) = b : 


x yy la? 


(5.31) 


r h , pVpM ] — p k ip il — p il p k i 

l"'x,x+li c x c x+l\ c x c x+l ^x'-'x+l- 


where b x {e) = —ies 3 G J. Identification of (5.24) with the Sutherland relation (5.30) is 
straightforward, based solely on the permutation action of Hamiltonian density 

(5.32) 

Multiplying the Sutherland relation by a string L x • • • L x _j from the left and a string 
L ;c+ 2 • • • L n from the right, summing over x and taking vacuum expectation value yields 
the defining relation for the amplitude operator 

[H, fl n ] = -ic (s 3 <8) O n _i - O n _i (8) s 3 ) , (5.33) 

where s 3 = e 11 — e 33 . Consequently, by expanding the unitary part of Liouvillian Cq, 


4(Poc) = i[H, Poo ] = i[H,Q n ]Sl - i n n [H,Q n }\ 


(5.34) 


in conjunction with (5.33), and employing the definition (5.20), the steady state 
condition (5.14) yields a decoupled system of boundary equations 

((vac| [v Al (U) - i(®S 1} - Bp)) = 0, 

(5.35) 


T>A 2 (B n ) + i 


| vac)) = 0, 


where two-leg boundary operators Bp, Bp € End (fH n 8> 77 a <8> 7f a ), reading 

3 3 

h e ij 


bx&'i ® la ® L j \ BP = ^2 ef b x <8> L ij <8> l a , 


(5.36) 


i,j =1 


i,j =1 


have been defined. Note that, due to (5.31), b x = i es x = —b x for £ e 


rat pairs of auxiliary operators (t + ,t ) and 


The last two lines of (5.24) indicate t 

(u + ,u“) span the Weyl-Heisenberg algebra. In conjunction with the highest weight 
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conditions (5.25) this fixes the representation of (t + ,t ) and (u + ,u ) to be that 


of a Fock space of two oscillator modes (bosons), specified by creation/annihilation 
operators, [bo-jb^.,] = 5 ay , [b^bo- /] = 0, a, a' G {t, |}, suggesting that the auxiliary 
space "H a is perhaps just a two-mode boson Fock space. While realization for all the 
other generators consistent with the bulk algebra g is not difficult to construct (e.g. v , 
P + F can be just the Schwinger boson representation of SU 2 - see 5th line of (5.24)), 


it turns out not to be consistent with the boundary conditions (5.25(.Therefore the 


auxiliary space "H a has to contain (at least) one additional degree of freedom. 


Ultimately, in order to fulfil (5.35), a straightforward calculation shows that it is 


enough to add a Verma module V p of complex spin representation (5.28) of s ( 2 and 
consider a triple-product space "H a = B ® B ®V P = lsp{ |j, k, l ); j, k, l G Z + }, and find a 


representation of the algebra (5.24) which is compliant with conditions 

/ | vac) 0 0 

L |vac) = I 7 /11,0,0) |vac) 0 

V 7 ?( 11 , 1 , 0 ) — | 0 , 0 , 1 » + 2 | 0 , 0 , 1 ) | 0 , 1 , 0 ) |vac) 


(5.37) 


(( vac 

( 1 , 0 , o| 

77((1,1,0| + ( 0 , 0 , 1 |)\ 


L = 0 

(vac 

^ 7 ( 0 , 1 , 0 | 

(5.38) 

V 0 

0 

(vac / 



with vacuum being given by the ground state |vac) = 10, 0, 0). These requirements 
are all satisfied by choosing representation (5.26|5.27|5.28) with p being fixed (5.29) as 


required by the conditions in the first two lines of (5.25). The last two lines of (5.25) 


are satisfied due to highest-weight-property of |vac). As such a representation is clearly 
irreducible, this concludes the proof . □ 


5.2. Grand canonical NESS and observables 


The formulae ( 5.16|5.19|5.23 5.29) yield explicit construction of a many-body density 
matrix of a family of degenerate NESSes = V^p 00 for any number of holes 
v G {0,1... n}. The computational complexity of obtaining any local information about 
the state poo, say to compute its matrix elements of the type (ii, ■ ■ ■ ,in\ Poo\ji, ■ ■ ■ ,jn) 
or local observables, is at most polynomial in n. Since the eigenspaces 'HN' of number- 
of-holes operator N 0 are orthogonal, one can also split decompose the Cholesky factors 

niT’OO = £ M n»00 

pM(s) = n<“ ) ( £ )S!<_■%), (5.39) 

since = o if 

v 7 ^ z/. Projected Cholesky factor satisfies a projected defining 


relation (5.33) 


[H, n^>] —ie (s 3 ®^,-^®- 

and can be expressed in terms of a constrained or microcanonical MPA 


QW = 


E « 


r \ (vac I L* 15 " 1 ■ • • L tnJn I vac) e 

[T, x Oi x ,2 ' 1 1 1 


nn 


olrijn 


(5.40) 


(5.41) 


^1 ijl •••In ijn 
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Note that since [f^, Nq] = 0, the Kronecker-5 constraint can just as well be replaced 


by 8 


(£, 


°ix, 2 


as only operators e lin 


o^njr 


for which Ex 2 = Ex Ey appear 


in MPA (5.19). 


We note two limiting cases of our NESS solution. For zero hole sector u — 0 one 
obtains exactly the fully polarized boundary driven isotropic (XXX) Heisenberg spin- 
1/2 chain and reproduces the solution reported in subsect. 2.5 The other extreme 
case (u = n) is the so-called dark state , i.e. a pure state /E _n ' ) = (e 22 )® n = 
12, 2 ... 2) (2, 2,... 21 which is unaffected by the dissipation, i.e. it simultaneously 
annihilated by Co and V, Copfo 1 = T>p^ = 0. 

Any convex mixture of states p^ = E^C/Po °v ^ is a valid NESS density 
operator as well, which factorizes (5.16) with a Cholesky factor Ct n = E*, x/c/EE 


Microcanonical constraint in (5.41) seems cumbersome as it prevents facilitating transfer 

There seems to be a particularly 
Namely, one may define a grand 


matrices for computation of local observables, 
attractive option which overcomes this problem. 
canonical nonequilbrium steady state (gcNESS) ensemble by taking a hole chemical 
potential p with c„ = exp (pu): 


Poo(e.m) = ^exp(fiv)p ( p(e). 


(5.42) 


u=0 


Note that the grand state corresponds to gcNESS with zero chemical potential poo(s) = 
Poo(e,p = 0). Clearly, the addition theorem for exponential function erases the 
constraint in MPO expansions: 


AO 

= . £ 

(vac L lljl (e, p) ■ 

■ • L lnjn (e,p) vac) e* Ul <g) • • • <g) e tnjn , 

(5.43) 






Poo (^) p) 

= £ 

((vac L tUl (e, p) ■ ■ 

■ ■ L* nJ ”(e, p)\v&c))e lljl <g) • • • <g) e lnJ ", 

(5,44) 


^1 ij 1 ■■•i'nijn 


where the chemical potential only modifies the components of the Lax operators as 

L l \e,p) = exp ( 7 ^, 2 ) L y (e), L lJ (e,p) = exp ( 7^,2 + ( 5 - 45 ) 

Moreover, introducing a transfer operator 

T(e,p) = ^L”(e,/r) = ^ L*(e, p) ® V j (e, p), (5.46) 

i i,j 

we define the grand canonical nonequilibrium partition function and express it via the 
transfer matrix method 


2n(E,p) = tr(poo(£,Ai)) = ({vac| (T(£,ft))’‘ |vac)). 


(5.47) 


The hole chemical potential p can be connected to the ensemble averaged filling factor 
(doping) r via logarithmic derivative of the partition function 


r := 


<E 


n 


ogZn(£ifl) . 

n Mv=o exp(n^t)trpoo 


(5.48) 


As usual, we expect the fluctuations (Sr) 2 = (-) 2 — r 2 to be thermodynamically small. 
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Expectation values of general (local) observables can again be extracted by 
facilitating the auxiliary vertex operators. Let Iuj = 0X0 lf^ n ^ be a 

generic local observable supported on a sublattice between sites x and y. Then, a 
formal expression 

(-^[x,j/]) Z n (O A 4 ) f' r ('^"[x,j/]Poo(^) aO)j (5.49) 

can be calculated from the MPA of p (e, /r) by tracing out the physical space TL n and 
associating to each observable X[ xy j a corresponding vertex operator via a mapping 
A £ : End (nf) -)• End('H a <g> ?Z a ), where £ = y — x + 1, using the prescription 

Ae(X) = X : = V tr (( e hjl 0 ■ ■ ■ 0 e idt )X ) V ljl ■ ■ ■ h ide . (5.50) 




For a complementary part of the lattice, i.e. where X[ X)V ] operates trivially, one has the 


transfer vertex operator T = Ai( 1 3 ), eq. (5.46), so the final expectation value reads 

(X [xM ) = Z~\{w aclT- 1 X T^lvac)). (5.51) 

For example, for on-site observables we have auxiliary vertex operators Ai(e* J ) = L :,z , 
e.g. for magnetization density Ai(s 3 ) = L 11 — L 33 . 

As for two point observables, we consider an interesting example of the current 
density tensor 

A 2 (J ij ) = ]S ij = i - h ij h ji ) (5.52) 


Stationarity (time-independence) of NESS and continuity equation (5.8) imply spatial- 


independence of current expectation values. In auxiliary transfer matrix formulation 


(5.46) this implies commutation of transfer vertex operator with current vertex operators 


when projected onto the subspace of states created upon action of T on the vacua, 
namely 

Ml[T,.FW>> = 0, Ml := «vac|T fc , |^» := T fc |vac». (5.53) 

Additionally, using the representation given in Theorem 2 and highest weight nature 
of the vacuum state, one can with some effort express the expectation values of total 


current operators (5.9) in terms of the nonequilibrium partition function (5.47) 

Zn-l 


(J 1 ) = 2e- 


M = 0, (J } = -2e 


Z, 


n—1 


(5.54) 


7, ’ ' ' z 

It would be challenging to attempt an analytic asmptotic computation of the grand 
canonical nonequilibrium partition function Z n {e,p) along the lines described in 

On the other hand, 
again, the 

_9 n n / • \ i n 

~ n 


subsection |3.2[ however this has not been accomplished yet. 
numerical computation c 
universal scaling Z n _ \ / Z n 


numerical computation of the expression (5.47) strongly suggests 

’~ 2 for all (generic) values of e, p. 
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6. Quasilocal conservation laws and linear response physics 

In this section we discuss the main ‘spin-off’ application of the exact solutions of 
boundary driven nonequilibrium master equations discussed so far, namely deriving 
novel, so-called quasilocal conserved quantities, and consequently deepening our 
understanding of the linear response physics of the corresponding closed (coherent, non- 
dissipative) models. The main exposition is focusing on the paradigmatic XXZ model, 
following Refs. [23176], while some comments with regard to other integrable chains will 
be given at the end. 


6.1. Universal R-matrix and exterior integrability of NESS density operator 


Using the concept of a universal R-matrix (see e.g. Refs. PI ESI 021 El ESI El mg) , one 
may find and explicitly construct the solution of YBE with U q (s\ 2 ) symmetry over an 


arbitrary triple tensor product of highest-weight Vernra modules (5.28) V. S| 

u Sl,S3 ( ( / J )-R'S2,S3 


V, 2 < 8 > V, 




R< 


■ 81 , 8 2 (<P - ^)Rsi,S 3 0)Rs2,S3 W = R S 2, S3 W R S 1, S3 (^) R S 1, S 2(V “ $) (6-1) 

for arbitrary representation parameters Si,s 2 ,S 3 G C and additive spectral parameters 
<p, d G C, where R S; acts nontrivially on the i-th and j-th module of the triple. For 
example Ri 1 (ip) = PR(^) is the six-vertex R-matrix of Eq. (2.39) up to permutation 


P of the pair of auxiliary spaces, while Ri s (ip) = L (<p,s) is a generic non-compact 


(non-Hermitian) Lax operator (2.40) so that YBE (6.1) for Si = S 2 — S 3 = s becomes 


the RLE relation (2.38). 


Taking the first two spaces as auxiliary and the third one as physical, V s <S)Vt( 8 )Vi = 
Ra ® Rb <8> R P , i.e., Si = s, s 2 = t, S 3 = |, and writing an infinite-dimensional, so-called 


exterior R-matrix as R a ,b(<£jS, t) = P a ,bR s ,i(^), Eq. (6.1) yields an alternative RLL 
relation swapping spectral and representation parameters in auxiliary tensor product of 
Lax operators 

Ra,b(v^ ^)La ,xilP 1 s)L h>x (0, t) L a a ,(d, t)Lb ja; (<^ 1 , s)R a b( < /? d, s, £). (6.2) 

Moreover, due to U q (s 1 2 ) invariance, the tensor product of highest-weight states, 
spanning a scalar representation, should always be left and right invariant under the 
R-matrix 

s. t) 10), |0)„ = |0>, |0) b , <0|. <0| b R,,b(v.», t) = <01. <0| b . (6.3) 

Remarkably, the relation (6.2) together with Q6.3 ) immediately implies commu¬ 
tativity of a two-parametric, highest-weight non-Hermitian transfer operator (HNTO) 
Wn(<p,s)e End(Rf") 

W n (<p,s) := (0| L(</p, s)® n |0), \W n (<p,s),W n {#,t)] = 0, (6.4) 

namely 
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= (°la (°lb |^n *) L b,*(</>, s )j R a,b(^ - $, S, t) |0) a |0) b = W n (d, t)W n (<p , s). 

Such a transfer matrix is essentially just the amplitude operator of NESS for 


a general asymmetric driving (2.91), specifically Q n (<p,s,x) = (</?, s)K(y/xY 


where the second, diagonal factor is inessential since it commutes with HNTO, 

W((M,*Vx)] = o. 

One may thus define m an exterior integrability of a nonequilibrium many-body 
density operator R^ = Ll n (<p, s, x)[ff n (<^, s, x)]t if Cholesky factor Ll n (ip,s,x) forms 
a commuting family for any values of driving parameters. Note, however, that as a 
manifestation of non-normality of the transfer operator W n (ip, s), \W n {p> 1 s ), Wff ($, t)] 7 ^ 
0 in general, hence the density operator R^p, s, x) does not form a commuting family. 

6.2. Non-Hermitian transfer operators with broken spin reversal symmetry and 
quasilocal conservation laws 


ffNTO (6.4) is neither a local operator, nor it is conserved in time as its time derivative 


is a non-local object, namely using ( 2.58[2.60 ) we End: 

[H, W n (ip, s)] = 2 sin q((cr z sin p sin r/s — < 7 ° cos p cos r/s) <E> W„_i (p, s ) 

— W n _i(tp, s) ® (cr z sin p sin r/s — a 0 cos <p cos r/s)). (6.5) 

Yet, it can be used to generate a very interesting family of operators in terms of 
differentiation with respect to the spin representation parameter s around the scalar 
point s = 0 

1 


Z n (p) = 


(sin p)' 


-d s W n (p,s)\ s=0 - {jj cot p)M*, 


( 6 . 6 ) 


where Mf = Xli—o ^ ® <7 Z <E> I 2 "- 1 -* is the conserved z—component of magnetization. 
The s —derivative can be implemented as MPA in terms of an additional ‘derivative 
anzilla’ qubit R c = C 2 , 

z.M = (°L {"WM 4 M • ■ ■ KM 10). |x>. - foot »>)*?. 

defining an extended Lax operator L '{p) e End(7f a <g) R c ® R p ) 

1 fL(p,0) <9 S L(<^, s)| s=0 ' 

L(<p, 0 ) 


L '{tp) = 


(6.7) 


= L 0 (^)l c + Li(p)af, ( 6 . 8 ) 


sm p \ 0 

where L 0 ((^) := (esc p)~L(p, 0), Li(<^) := (esc <^)<9 S L((^, s)| s=0 . We shall refer to the 
operator family Z n (p) as the modified highest-weight non-Hermitian transfer operators 


(mffNTO). Note that the Z-operator (2.113) yielding the first-order perturbative 
expression of NESS is just Z = Z n (7 t/2 ) t . It can be shown [76] that in the massless 
regime (for a dense set of real 17 ) Z n (tp ) are quasilocal operators whose time-derivative 
is localized at the chain boundaries for a suitable domain p e T> C C. Indeed, 


differentiating (6.5) w.r.t. s at s = 0 and using the definition ( 6 . 6 ) we immediately 
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obtain a very insightful relation 

[H, Z n (<p)\ = 2/7 sin 77 (<r z <g) t 2 n-i — t 2 n-i <8> a z ) 

— 2 sin r] cot ip (cr° (8) Z n _\{ip) — Z n _\(ip) <8> cr°) . (6.9) 

Writing the Lax operator components L Q G End('H a < 8 > 'HQ, L" G End('H a ), via 
L '{ip) = L a (^) <8> cr“, L ,a (<£>) = Lj} (<^)1 c + L"(<^)cr+ satisfying boundary transition 

conditions 


<o|.<oicL' 0 = <oL<o| c , 
L'°| 0 )J 1 ) =| 0 )J 1 ) C , 


,(0| c L'+=0, 
|0>.|l>o = 0. 


L'“|0).|l)„ = »cot V |0)J0) c , L z,± |0) a |0) c — 0, (6.10) 

one sees that mHNTOs allow for an expression in terms of open boundary translationally 
invariant sum of local operators 

n n—r 

Z n(<f) = !fl_2 x ^ Qr(^p') ^ H 2 n— r —Xj (6.11) 

r =2 x=0 

in terms of local r— point densities q r {<p) G End('Hp r ) with explicit MPA representation, 
which is obtained by careful inspection of the definitions ( 6.7[|6.8 ) 

?■■(¥>)= E (oi L o( 6 . 12 ) 

Q2...Q r -l£j’ 


Using the local operator sum ansatz (6.11) one is able to rewrite the RHS of (6.9) in a 
form of a sum of operators localized at the boundaries 

[H, Z n (<p)\ = 2?/sin // (<r z (8> l 2 «-i — t 2 n-i (8) cx z ) 

n 

+ 2 sin rj cot (p (g r (</?) <E> l 2 n ~ r ~ <S> g r (<^)) • (6.13) 

r=2 

We shall now demonstrate (following [761 [75] ) that there are important parameter 
regimes for which the operator sequence {q r { ( p)]T' = 2,3...} is quickly decreasing in 
a suitable operator norm, so the operator family ( 6 . 11 ) can be considered as quasilocal 
and almost conserved. Moreover, the operators Z n (ip) are represented as strictly 
lower triangular matrices with zero diagonal, (v\ Z n (ip) \vf) = 0 if ord(z/) < ord(z/), 
(u\ Z n (ip) | u) = 0 for any ip. Hence they are non-diaginalizable for any n > 1 as their 
spectrum contains only 0 . 

However, in order to formulate our results precisely we first state two definitions of 
‘essential’ locality of translationally invariant operators: 

Quasilocality: An operator sequence Z n G End("Hp n ) which can be written as an open 
boundary translationally invariant sum of local operators q r , like ( 6 . 11 ), for any n, is 
called quasilocal if there exist positive constants 7 , £ > 0 , such that 

IMIhs < (6.14) 

where, for any matrix a, 

2 tr(eda) 


IHS ■- 


trl 


(6.15) 
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is a normalized Hilbert-Schmidt norm which satisfies a nice extensivity property 

IMIhs = || a <£> 1 cz||hs, Vd, (6.16) 

as well as the normalized Cauchy-Schwartz inequality 
tr(a&) 


trl 


< H hs 


|HS- 


(6.17) 


We remark that the Hilbert-Schmidt operator norm is the natural norm for high- 
temperature statistical mechanics as it is linked to an infinite temperature, tracial state 
cuo(a) = tra/trl, specificallly 11a11 = <^o {a^a). Note also that it satisfies a useful 
inequality in relation to a C* operator norm ||&|| 2 = sup^cnamely for any pair of 
bounded operators a,b (say, elements of End(7fp n )), ||a6||ns < ||ci||hs||&|I- 
Pseudolocality: An operator sequence Z n 6 End('Hp n ) of the form (6.11) is called 
pseudolocal if there exists a positive constant K > 0, such that 

IIHhs — An. (6.18) 

Clearly, quasilocality implies pseudolocality as follows straightforwardly from the 
definitions. In order to demonstrate locality of rnHNTO Z n (p) for XXZ one needs 
to study the sequence of Hilbert Schmitd norms ||^(v 7 )IIhs- In fact; f° r a set of 
commensurate anisotropies r] = —, l,m € Z + , densely covering the easy-plane regime 
|A| < 1, one can explicitly study a general inner product US 

K r(<p,<p') ■-= ftr(qj((p)q r (</))=( —;) (1 T(ip, ip') r ~ 2 |1), r > 2. (6.19) 

in terms of iterating the following finite, m —dimensional transfer matrix acting over a 
vector space lsp{|/c), k — 1,..., m}, 

771—1 m— 2 


T{p,p')=>2(cl + COtpCOtp'sl) l fc )( fc l + Y1 


| 'Sfc'Sfc+l | 


7 (\k){k+l\ + l^ + l)(^l); (6.20) 


2 sin p sin p 1 

k= 1 fc=l r r 

where c*, := cos(7r lk/m),Sk ■= sin(7r Ik/m). Straightforward calculation [76, 75j shows 
that the leading eigenvalue r = e~ 2 ^ of T(tp,tp') has modulus smaller than 1, i.e. 
£ > 0, and hence Z n (p) being quasilocal, exactly if p, p' belong to the vertical strip 
V m = { p ; |R eV 9 — || < Moreover, on V m one can easily compute the full extensive 
normalized Hilbert-Schmidt inner-product 
, , tiA*B 

( ; ) - trl ’ ( 6 ‘ 21 ) 

which turns the space of observables End(H® n ) into a Hilbert space, between all 
rnHNTO. Namely, for any pair p, p’ G T> m 

n 

(Z n (p),Z n (p')) = J^(n-r + l )n r (p,p') 


r=2 
00 


= n 


y K r (p, p') - ^(r - 1 )n T {p, p') + 0(ne 


-2£n\ 


r=2 


r=2 
0\ 


= nK(p,p') + O(7i 0 ),, 
{Z T n {p),Z n {p')) =0, 


( 6 . 22 ) 
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where 

K M) = X> r (</h¥>')= ( sin ^ iy/ ) (1| (1 — T(^>, |1) 

8 V 2 sin((m — l)(y? + <p')) 
sin ip sin <// sin (m(<p + ip 1 )) 

With a simple technical trick Eza l6lj one can, again for commensurate anisotropies 
7] = nl/m densely covering the gapless regime |A| < 1, define a set of quasilocal 
conserved operators which are exactly conserved for the XXZ Hamiltonian with periodic 
(or even twisted ESI) boundary conditions, 

Hpbc — H+2a + ®t- 2 n-2®cr~ +2a~ ®t- 2 n-2®a + +Aa z Z>t2n-2®a z , (6.24) 

namely 

Yn(v) = , . 1 s n dsK(<l,s)|s =0 - (r) cot tp)M Z 
(sm (p) n 

= tr a (0| c Li(^)L' 2 (</?) • ■ • L ’ n (<p) |l) c - (rj cot ip)M z , (6.25) 

where 


V n (ip, s ) = tr a {Li(</?, s)L 2 (ip, s) ■ ■ • L„(</?, s)} , (6.26) 

with the auxiliary space being naturally truncated to "H a = lsp{|/c), k = 0,1..., m — 1} 
since (m — 1|L _ (^, s) = 0. The Sutherland condition (2.47) then immediately implies 
[Hpbc, Vn(<p, s)] = 0, and consequently [via (6.25)], exact conservation 


[H phc ,Y n {(p)\ = 0, VV e C, 


while YBE (6.2) implies 


\Y n {ip),Y n (ip')\ = 0, VV, <p' G C. 


(6.27) 

(6.28) 


Furthermore, one can easily show m that the difference between rnHNTO Z n (ip) 
and the so-called modified periodic non-Hermitian transfer operator (rnPNTO) Y n (p) 
is exponentially small away from the boundaries. More precisely, rnPNTO is again 
quasilocal in the sense that it can be written as a translationally invariant sum of local 
operators (with the same densities as Z n (ip)) 

n n— 1 

Y n(x) = <S x (l 2 n-i r 0 q r (<p)) + y n (<p) (6.29) 

r =2 *=0 

where S : End("Hp") —> End('Hp n ) is a left-shift rotation map which is completely 
specified by the action on the Pauli basis 

S(a Q ° 0 cr ai 0 ■ ■ ■ a“™- 2 0 a"- 1 ) = a ai 0 a" 2 0 • • • a"- 1 0 a“°, (6.30) 


and the remainder is exponentially small in Hilbert-Schmidt norm ||2/n(<^) ||hs = 0(e~Z n ). 
Similarly to HNTO (6.22), the family of rnPNTO has asymptotically the same kernel 
of inner products 


(U(0), W)) =»*W) + 0(e- 2 «"), {Ypt),Y„(v')) = 0(e (6.31) 
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The standard algebraic Bethe ansatz (ABA) machinery [45[ 251 [29] allows one 
to derive a sequence of strictly local translationally invariant and exactly conserved 
operators Q r , r = 2 , 3 .. [Q r , Q r /\ = 0, 

n— 1 


Qr=dp' tog K(V,§)I».= 3 = 5^(1 


^n — r 


Q 


(r 




(6.32) 


x=0 


with the first term in the series being proportional to the hamiltonian Q 2 oc H p b c , 
and where q^ r ’ G End('Hp r ) are the corresponding local densities. Importantly, the 


ABA transfer operator V n (ip, |) is spin reversal (2.19) invariant, and hence are all local 
conserved operators 

PV n {<p, DP - 1 = V n {<p, §), PQrP - 1 = Qr. (6.33) 

On the other hand, the non-Hermitian transfer operators, and the corresponding 
quasilocal conserved operators, satisfy a more specific PT-likc symmetry (analogous 


to the one discussed in [72]). following from (2.52) 

PZ n (r)P-'=Zl (ir-p), PY n (v)P- l =Yp x-ip). (6.34) 

This means that the modified (and quasilocal for 77 = 7 rl/m) transfer operators can be 
decomposed into even and odd w.r.t. spin reversal 

Ztw = Z»M ± PZn(v)P~\ Y*(v) = ± PY n ( V )P~\ (6.35) 

such that 

[i/ pbc , Y‘(r)\ = 0, PYMP- 1 = (6.36) 

«,“(■?), U"'(»5)) =nS a , a .K( v , v ’)+0(e-^\ (6.37) 

а, a' e {±}, and similar relations for Zf((p) with open boundaries. 

б. 3. Lower bounds on high temperature ballistic transport coefficients 

Existence of quasilocal conserved operators is extremely interesting for deriving bounds 
on linear-response transport coefficients. For example, considering the extensive current 

71—1 

J = ^dP(l 2 »-2 <g> j) (6.38) 

x=0 

(with periodic boundary conditions), the famous Green-Kubo formula expresses the 
corresponding conductivity a' (uj) (in our case spin-conductivity if j is a local spin 


current (3.6), but for a general discussion J can be any current linked to an appropriate 


conservation law) in terms of the current auto-correlation function. In particular, 


= lim lim — [ dt'e lut (J(t'), J{0))g 
t—¥ oo n—¥ oo Tl J q 


(6.39) 


where J{t) := e lHphct Je lHpbct is the Heisenberg dynamics of the corresponding current 
operator, and 


1 rfl +-r» ( P ^pbc Rg (/^ ^)^pbc^ 

AX tre-W*. - L 


(6.40) 
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is the Kubo-Mori inner product, which reduces to Hilber-Schmidt inner product (6.21) 
in the limit of infinite temperature, lim^c^A B)p = (A,B). When d.c. (c o = 0) 
conductivity diverges, one defines the (spin) Drude weight D 

cr'(u) = 2ttD5(uj) + cr reg (a;) 

which can be again expressed with a Green-Kubo-like formula 


(6.41) 


D= lirn lim [ d t'(J(t'), J(0))p. 

t ^oo n—>oo Ztn /n 


(6.42) 


Drude weight can be considered as a ballistic transport coefficient, and D > 0 signals 
an ideal, ballistic spin transport in an extended system at finite temperatures 
At high temperature fd —> 0, the leading order Drude weight can be expressed as 


D = (3D OC + O(l3 2 ), D 0 


1 

= lim lim -— / dt'(J(t'), J). 

t-> oo n^-oc Ztn In 


(6.43) 


For integrable quantum systems, having a (countable) set of local extensive conserved 


operators {Q r } (6.32), Zotos, Naef and Prelovsek [ 100] suggested to use Mazur bound 
ra. to rigorously estimate the high-temperature Drude weight from below, 


Doc > lim — y^(J,Q r )(K 1 ) r y(Qr',J) 

n —vno Zn ZJ 


(6.44) 


rr' 


where K r y := ( Q r ,Q r >) is a positive-definite matrix of inner-products of independent 
conserved operators. Usually, one picks such linear combinations of Q r which are 
mutually orthogonal, say, using a Gram-Sclnnidt procedure, so the above formula 
simplifies with (K~ 1 ) r y = 5 r y ||<5rIlni- 

Howe ver, the situation becomes interesting and quite intricate for systems with 
Z 2 symmetries, like spin-reversal, parity-hole, etc, such that the corresponding spin or 
charge current is odd under the transformation 

PJP- 1 = -J, (6.45) 

and the symmetry in the corresponding equilibrium state is un-broken, tr(Xe~ /3H ) = 
tr (PXP- 1 e~P H ),VX, i.e. in the absence of external magnetic fields, chemical potentials, 
etc, such that PHP -1 = H. In the case of XXZ model this immediately implies that 
the RHS of Mazur bound (6.44) has to vanish, since ( J,Q m ) = —(-P JP 1-1 , PQmP 1 ) = 
— (</, Qm) — 0, and one has to rely on effective theories and approximations [S, 85] . 


However, the situation drastically changes due the presence of quasilocal conserved 
operators with odd spin reversal symmetry {Y~ (</?)}■ Even replacing a single operator 
from this set Y~(n/2) for commensurate anisotropy 77 = nl/m (or, equivalently, Z n (<p/2) 
for open boundaries [3S]) into RHS of (6.44) one obtains a non-vanishing lower bound 
[69] [35] [6l] 


Doc > D z = sin 2 (7rZ/m) 


m 


m — 1 


(6.46) 
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Figure 9. (from [76]) Optimized Mazur lower bound on high-temperature spin Drude 
weight Ck (black) [75] versus the less optimal bound C\ of Ref. [55] (red) which is 
based on a single quasilocal almost-conserved operator Z n (ip = 7r/2) (or Y n ( 7r/2)), for 
XXZ model as a function of anisotropy parameter A. 


However, one can do much better than that! Considering the full continuous family of 
non-Hermitian quasilocal conserved operators {Y n (cp),cp G V m } replacing a countable 


family { Q m , m — 1, 2 ...} in (6.44) on can write the Mazur bound 


D 


oo > D 


K A 


> D z , D k = -Re / d 2 <pa{<p)f(<p), 


(6.47) 


’Vr, 


where a{ip) := (j XiX +i,Y n (<p)), in terms of a solution /(</?) of a complex Fredholm 


equation of the first kind, involving the kernel (6.23) (see Ref. [75] for details) 
1 
2 




(6.48) 


I T>„ 


In our case a(ip) = i/4, and the Fredholm equation has a simple explicit solution 
/(</?) = —-fms'l\ sin</?|~ 4 , yielding an explicit expression for the optimised Mazur bound 

' 2 (7t l/m) m . (2tt\ 


D k — 


sm 


sin 


(n/m) 


1-sin — 

2n \ m J 


(6.49) 


It could be tempting to speculate that this bound is in fact saturating the high- 
temperature spin-Drude weight (see Fig. [9]), in particular since it agrees very well with 
the state-of-the-art time dependent DMRG simulations B5- Note, however, that both, 
the non-optimal and optimised lower bounds Dz, Dk are no-where continuous (‘fractal’) 
functions of the anisotropy parameter A = cos 7]. This suggests an interesting option 
that thermodynamic properties such as transport coefficients of strongly interacting 
systems might be no-where continuous functions for a finite range of parameters. 
Similar analysis could be attempted for finite inverse temperatures /3 > 0, but then the 
evaluation of the K(ip,ip') and the function a(ip) should be evaluated either numerically 
or approximately (perturbatively). 
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6.f. Models with non-deformed symmetries and lower bounds on high temperature 
diffusion constants 


Analysis of the previous section showed how new quasilocal operators can be distilled 
from NESS density operator in boundary dissipatively driven model with trigonometric 
A-matrix (corresponding to U q {s\ 2 ) quantum symmetry with uni-modular q). We expect 
that similar results could be obtained for other integrable models sharing the same A- 
matrix, say Sine-Gordon held theory in 1 + 1 dimensions, or its discretisation, the 
quantum Hirota equation (see e.g. 12a ). However, our conserved operators Y n (cp) or 
Z n (tp) are no-longer quasilocal when \q\ > 1, when the amplitudes of the MPA diverge 
super-exponentially, or even in the un-deformed slo symmetric limit q —> 1. 

In such a case, for the XXX spin 1/2 model, the central quasilocal operator 
Z T (ip = 7t/ 2) goes to 

n— 1 n 

Z = -i<9 £ fi n | £=0 = a x a y’ ( 6 -50) 

x=l y=x -\-1 

so it becomes quadratically extensive in the sense that ||Z||^ S — \qn 2 with q being some 
real constant, while still satisfying almost conservation, or conservation law property 
[H, Z] = — o\ + cr^. One hnds very similar behaviour for other models with Lie 
(underformed) symmetries presented in this overview, namely Fermi-Hubbard and spin- 
1 Lai-Sutherland models. 

For example, for symmetrically boundary driven Hubbard model, one can define 
an analogous operator as 

Z = —i<9 £ fl n (A = 0, w = |e) = (6.51) 

n—l x<y / y -1 \ / y-1 \ 

+ t x t x+i) - 2 u ' 52 (- i ) x ~ v<7 t (n) v t * (n r " ) r y~’ 

which is again quadratically superextensive ||Z||g S ~ \qn 2 and satishes almost 
conservation property 


[H, Z] = + < + rl 


(6.52) 


Note that both operators, (6.50) and (6.51), can be viewed as level-1 generators of the 
Yangian symmetry of the respective models, truncated to a finite size n 0193]. Similarly 
as in the case of XXZ model, the entire two-parameter set of operators H n (A, w) (4.39) 
can be considered as an HNTO, namely one hnds by explicit computation that 

[O n (A, w), f2 n (A', it/)] = 0, VA, G C. (6.53) 


However, proving the existence or even explicitly constructing the corresponding 
intertwiner, the (infinite-dimensional) exterior A-matrix, remains a problem for the 
future. 

In the absence of local and quasilocal operators which are odd under spin reversal, or 
particle-hole transformation [ffl P , one may attempt use the almost-conserved Hermitian 


ff Reader should remember that the ABA transfer operator and all the derived local conserved operators 
are even under P. 
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operator Q = \(Z — Z t), PQP _1 = — Q to estimate d.c. spin transport coefficients. 
However, due to super-extensivity, the contribution of Q to Drude weight is vanishing 
in TL. Nevertheless, using a careful estimation of the effects of time evolution of the 
derivative [H, Q] (which is by assumption of almost-conservation localised near the 
boundaries) by means of the Lieb-Robinson bound pTJ, one can find a rigorous lower 
bound on the Green-Kubo expression for the high-temperature diffusion constant m 



(6.54) 


For the locally interacting hamiltonian p[ = 'Yh x h X)X + 1 , local current observable 
J = J2 x jx,x +1 and general quadratically extensive almost conserved operator Q, with 
q = lim^oo ~ 2 11 Q 11 pjs > 0, one finds a general theorem [74], stating that 



(6.55) 


where v is the Lieb-Robinson velocity [32., 58], which can be estimated in terms of the 
norm of local Hamiltonian density, v < 6||/i||, hence D diS > \(j x , x + 1 , Q)\ 2 /mh\\q)- 
Applying to XXX and Fermi-Hubbard models this result implies strict bounds on the 
respective spin and charge diffusion constants, D x ^ x >1/6, H^ i g bbard > l/(3w 2 ), which, 
from a rigorous point of view, prove that the high-temperature spin/charge transport 
in these models cannot be sub-diffusive or even insulating. Whereas in the Heisenberg 
model this bound may be superfluous, as DMRG numerical simulations suggest that the 
high-temperature spin transport in the isotropic point seems to be anomalous (super- 
diffusive but sub-ballistic, and hence D diS = oo) [98] . in the Hubbard model the bound 
seem to be less trivial as the numerics suggests diffusive transport Hi. 

7. Discussion 

7.1. Open problems 

We shall close the presentation of this growing subject with a list of, to author’s taste, 
most urgent open problems. 

• We see currently no analytical technique to compute the Liouvillian spectrum, its 
gap and decay modes, i.e. to solve the full Liouvillian eigenproblem Cvj = XjVj in 
the models where NESS (fixed point) is an exactly solvable MPA. Even in the 
simplest lion-trivial (strongly-interacting) case of integrable NESS, e.g. in the 
boundary driven XXX spin 1/2 chain, we currently do not understand how to build 
higher decay modes Vj, with decay rates Xj, Re Xj < 0 (see e.g. Refs. [67], [68, [54] 
for such results on non-interacting systems). It is not even clear at present if the 
integrability of NESS implies that the problem of diagonalizing the full Liouvillian 
C needs to be integrable. One should perhaps note that non-trivial statements 
can be made about the structure of decay mode spectrum based on rather general 
properties of the Liouvillian such as an analogue of the topical PT-symmetry ra- 
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• Alternatively, one may try to construct exact solutions for relaxation dynamics p{t) 
directly in the time domain for specific non-trivial initial states p(0), in analogy with 
SEP [82]. For example, one may devise a quench protocol, where p(0) is an exact 
MPA NESS of an integrable model. Then, at t — 0, one suddenly changes the 
parameter of the Hamiltonian or of the dissipator, such that the NESS of the after- 
quench problem remains integrable. It is perhaps reasonable to expect that then 
the full dynamics p(t ) remains integrable, i.e. exactly solvable, as well. 

• So far, one is able to write exact MPA solutions of NESS only for integrable quantum 
chains with very specific dissipative boundary conditions, which can be phrased as 
a pure source and a pure sink. More general boundary conditions can be treated 
only perturbively in the system-bath coupling constant. This situation is quite 
different than in SEP P, 82] where one can typically exactly treat the most general 
local boundary conditions. One perhaps needs to extend the Skylanin’s concept of 
reflection algebra [8.6] to quantum Liouvillc space formalism. 

• The NESS density operator p^ with exact MPA structure could be compared with 
its equilibrium integrable counterpart, which is the Gibbs operator e~^ H where H 
is a Hamiltonian of an integrable quantum chain, logpoo can thus be considered 
as a kind of nonequilibrium integrable Hamiltonian and the eigenvalues of can 
be considered as probabilities. More generally, p^ determines the nonequilibrium 
thermodynamic state of the system and one may be interested in computing its 
observable properties. Thus it would be desirable to have a Bethe-ansatz for 
diagonalization of p^. In Ref. 177] the first step of such protocol has been outlined, 
i.e. the single quasi-particle spectrum of p^ for boundary driven XXX spin 1/2 
chain has been calculated, but problems with higher quasi-particle excitations have 
been identified. 

• The NESS density operator p^ only entails average steady-state properties of the 
system. In order to access fluctuation properties, such as e.g. cumulants of the 
current, one needs to go beyond the master equation and consider the so-called 
full-counting-statistics ra or analogous large-deviation-theory formalism [92] . It 
has been shown in Ref. [3T|, that the fc-th cumulant problem is exactly solvable by 
MPA for classical ASEP for any k. Also, for a boundary driven quantum XXZ 
spin 1/2 chain we were able to compute all cumulants of the current perturbatively 
in the system-bath coupling strength. However, an analog quantum problem to 
Ref. [31] . i.e. computing exact current cummulants in boundary driven XXZ spin- 
1/2 chain, remains open. A partial, perturbative, result in this direction has been 
achieved in Ref. ra, where all cumulants of the current have been calculated in 
the lowest two orders of system-bath coupling. 

• All exact NESS solutions presented in this topical review refer to one-dimensional 
chains with ultra local dissipation acting only on the first and the last site of the 
chain. It would be tempting and physically desirable to extend our nonequilibrium 
integrability techniques to quantum held theories with incoherent particle sources 
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and sinks at or near the boundaries. Prime candidates are Sine-Gordon and Lieb- 
Liniger models, highly relevant for low-energy condensed matter or cold atom 
physics. The main difficulty is entailed in regularisation of the dissipator for a 
field theory, namely it should correspond to a source/sink of a particle with a well 
defined single-particle wave-function [S8] which cannot be located strictly at a point 
since this would correspond to infinite energy. 

• Motivated by numerical results of Ref. [80], showing a phase transition from ballistic 
to diffusive spin-transport in a classical lattice Landau-Lifshitz spin chain (see 
e.g. [26]) which can be considered as an integrable classical limit of the XXZ 
chain, one may be tempted to formulate a consistent classical limit for integrable 
nonequilibrium boundary driven models (i.e. ‘classical exterior integrability’). For 
example, one may write a boundary driven lattice Landau-Lifshitz model with 
Langevin noise processes attached to the end sites and attempt to construct an 
exact MPA for the steady state (classical NESS). In full analogy with the story 
on XXZ spin 1/2 chain, one again has a spin reversal symmetry, with respect to 
which all classical local conserved quantities are even, but one expects to derive new 
quasilocal conserved quantities with broken spin-reversal symmetry which could 
explain the ballistic classical spin transport. 

• We note that an alternative path to integrability in open nonequilibrium systems in 
terms of scattering state formalism, which is specially suited for quantum impurity 
problems , has been outlined in Ref. [55] • It is not clear whether and how a link 
to the boundary dissipation approach discussed in this review can be made, and 
whether the latter could be implemented to treat integrable impurity problems. 

1.2. Conclusion 

This topical review presented a state-of-the-art (or better to say, a snap-shot in 
developing the theory of) exact MPA solutions of steady states of dissipatively boundary 
driven quantum integrable chains. An attempt to make a coherent presentation covering 
a variety of different integrable models under the same footing has been made. The key 
constructive (algebraic) methods have been identified and related to general methods of 
quantum integrability, such as the Lax structure and Yang-Baxter equation. However, 
important distinction to integrable equilibrium problems should be underlined, namely 
in dissipatively driven quantum chains one should consider non-unitary irreducible 
representations of the quantum (deformed), or Lie symmetries of the model, where the 
representation parameter of these algebraic structures is connected to the noise strength 
at the chain ends. 

Apart from reviewing previously published material in a coherent and self-contained 
manner, this article contains also several original scientific results, the most notable 
being: (i) MPA for NESS in asymmetrically driven XXZ chain and with arbitrary 
(asymmetric) transverse fields at the boundary, and (ii) exact asymptotic computation 
of the nonequilibrium partition function for the XXX model (but also extending 
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to other models with un-deformed Lie symmetries) which is the basis for computing 
nonequilibrium thermodynamics and observables. 
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